LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
zstemr.f
Go to the documentation of this file.
1 *> \brief \b ZSTEMR
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZSTEMR + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zstemr.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zstemr.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zstemr.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
22 * M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK,
23 * IWORK, LIWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBZ, RANGE
27 * LOGICAL TRYRAC
28 * INTEGER IL, INFO, IU, LDZ, NZC, LIWORK, LWORK, M, N
29 * DOUBLE PRECISION VL, VU
30 * ..
31 * .. Array Arguments ..
32 * INTEGER ISUPPZ( * ), IWORK( * )
33 * DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
34 * COMPLEX*16 Z( LDZ, * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> ZSTEMR computes selected eigenvalues and, optionally, eigenvectors
44 *> of a real symmetric tridiagonal matrix T. Any such unreduced matrix has
45 *> a well defined set of pairwise different real eigenvalues, the corresponding
46 *> real eigenvectors are pairwise orthogonal.
47 *>
48 *> The spectrum may be computed either completely or partially by specifying
49 *> either an interval (VL,VU] or a range of indices IL:IU for the desired
50 *> eigenvalues.
51 *>
52 *> Depending on the number of desired eigenvalues, these are computed either
53 *> by bisection or the dqds algorithm. Numerically orthogonal eigenvectors are
54 *> computed by the use of various suitable L D L^T factorizations near clusters
55 *> of close eigenvalues (referred to as RRRs, Relatively Robust
56 *> Representations). An informal sketch of the algorithm follows.
57 *>
58 *> For each unreduced block (submatrix) of T,
59 *> (a) Compute T - sigma I = L D L^T, so that L and D
60 *> define all the wanted eigenvalues to high relative accuracy.
61 *> This means that small relative changes in the entries of D and L
62 *> cause only small relative changes in the eigenvalues and
63 *> eigenvectors. The standard (unfactored) representation of the
64 *> tridiagonal matrix T does not have this property in general.
65 *> (b) Compute the eigenvalues to suitable accuracy.
66 *> If the eigenvectors are desired, the algorithm attains full
67 *> accuracy of the computed eigenvalues only right before
68 *> the corresponding vectors have to be computed, see steps c) and d).
69 *> (c) For each cluster of close eigenvalues, select a new
70 *> shift close to the cluster, find a new factorization, and refine
71 *> the shifted eigenvalues to suitable accuracy.
72 *> (d) For each eigenvalue with a large enough relative separation compute
73 *> the corresponding eigenvector by forming a rank revealing twisted
74 *> factorization. Go back to (c) for any clusters that remain.
75 *>
76 *> For more details, see:
77 *> - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations
78 *> to compute orthogonal eigenvectors of symmetric tridiagonal matrices,"
79 *> Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
80 *> - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and
81 *> Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25,
82 *> 2004. Also LAPACK Working Note 154.
83 *> - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
84 *> tridiagonal eigenvalue/eigenvector problem",
85 *> Computer Science Division Technical Report No. UCB/CSD-97-971,
86 *> UC Berkeley, May 1997.
87 *>
88 *> Further Details
89 *> 1.ZSTEMR works only on machines which follow IEEE-754
90 *> floating-point standard in their handling of infinities and NaNs.
91 *> This permits the use of efficient inner loops avoiding a check for
92 *> zero divisors.
93 *>
94 *> 2. LAPACK routines can be used to reduce a complex Hermitean matrix to
95 *> real symmetric tridiagonal form.
96 *>
97 *> (Any complex Hermitean tridiagonal matrix has real values on its diagonal
98 *> and potentially complex numbers on its off-diagonals. By applying a
99 *> similarity transform with an appropriate diagonal matrix
100 *> diag(1,e^{i \phy_1}, ... , e^{i \phy_{n-1}}), the complex Hermitean
101 *> matrix can be transformed into a real symmetric matrix and complex
102 *> arithmetic can be entirely avoided.)
103 *>
104 *> While the eigenvectors of the real symmetric tridiagonal matrix are real,
105 *> the eigenvectors of original complex Hermitean matrix have complex entries
106 *> in general.
107 *> Since LAPACK drivers overwrite the matrix data with the eigenvectors,
108 *> ZSTEMR accepts complex workspace to facilitate interoperability
109 *> with ZUNMTR or ZUPMTR.
110 *> \endverbatim
111 *
112 * Arguments:
113 * ==========
114 *
115 *> \param[in] JOBZ
116 *> \verbatim
117 *> JOBZ is CHARACTER*1
118 *> = 'N': Compute eigenvalues only;
119 *> = 'V': Compute eigenvalues and eigenvectors.
120 *> \endverbatim
121 *>
122 *> \param[in] RANGE
123 *> \verbatim
124 *> RANGE is CHARACTER*1
125 *> = 'A': all eigenvalues will be found.
126 *> = 'V': all eigenvalues in the half-open interval (VL,VU]
127 *> will be found.
128 *> = 'I': the IL-th through IU-th eigenvalues will be found.
129 *> \endverbatim
130 *>
131 *> \param[in] N
132 *> \verbatim
133 *> N is INTEGER
134 *> The order of the matrix. N >= 0.
135 *> \endverbatim
136 *>
137 *> \param[in,out] D
138 *> \verbatim
139 *> D is DOUBLE PRECISION array, dimension (N)
140 *> On entry, the N diagonal elements of the tridiagonal matrix
141 *> T. On exit, D is overwritten.
142 *> \endverbatim
143 *>
144 *> \param[in,out] E
145 *> \verbatim
146 *> E is DOUBLE PRECISION array, dimension (N)
147 *> On entry, the (N-1) subdiagonal elements of the tridiagonal
148 *> matrix T in elements 1 to N-1 of E. E(N) need not be set on
149 *> input, but is used internally as workspace.
150 *> On exit, E is overwritten.
151 *> \endverbatim
152 *>
153 *> \param[in] VL
154 *> \verbatim
155 *> VL is DOUBLE PRECISION
156 *> \endverbatim
157 *>
158 *> \param[in] VU
159 *> \verbatim
160 *> VU is DOUBLE PRECISION
161 *>
162 *> If RANGE='V', the lower and upper bounds of the interval to
163 *> be searched for eigenvalues. VL < VU.
164 *> Not referenced if RANGE = 'A' or 'I'.
165 *> \endverbatim
166 *>
167 *> \param[in] IL
168 *> \verbatim
169 *> IL is INTEGER
170 *> \endverbatim
171 *>
172 *> \param[in] IU
173 *> \verbatim
174 *> IU is INTEGER
175 *>
176 *> If RANGE='I', the indices (in ascending order) of the
177 *> smallest and largest eigenvalues to be returned.
178 *> 1 <= IL <= IU <= N, if N > 0.
179 *> Not referenced if RANGE = 'A' or 'V'.
180 *> \endverbatim
181 *>
182 *> \param[out] M
183 *> \verbatim
184 *> M is INTEGER
185 *> The total number of eigenvalues found. 0 <= M <= N.
186 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
187 *> \endverbatim
188 *>
189 *> \param[out] W
190 *> \verbatim
191 *> W is DOUBLE PRECISION array, dimension (N)
192 *> The first M elements contain the selected eigenvalues in
193 *> ascending order.
194 *> \endverbatim
195 *>
196 *> \param[out] Z
197 *> \verbatim
198 *> Z is COMPLEX*16 array, dimension (LDZ, max(1,M) )
199 *> If JOBZ = 'V', and if INFO = 0, then the first M columns of Z
200 *> contain the orthonormal eigenvectors of the matrix T
201 *> corresponding to the selected eigenvalues, with the i-th
202 *> column of Z holding the eigenvector associated with W(i).
203 *> If JOBZ = 'N', then Z is not referenced.
204 *> Note: the user must ensure that at least max(1,M) columns are
205 *> supplied in the array Z; if RANGE = 'V', the exact value of M
206 *> is not known in advance and can be computed with a workspace
207 *> query by setting NZC = -1, see below.
208 *> \endverbatim
209 *>
210 *> \param[in] LDZ
211 *> \verbatim
212 *> LDZ is INTEGER
213 *> The leading dimension of the array Z. LDZ >= 1, and if
214 *> JOBZ = 'V', then LDZ >= max(1,N).
215 *> \endverbatim
216 *>
217 *> \param[in] NZC
218 *> \verbatim
219 *> NZC is INTEGER
220 *> The number of eigenvectors to be held in the array Z.
221 *> If RANGE = 'A', then NZC >= max(1,N).
222 *> If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU].
223 *> If RANGE = 'I', then NZC >= IU-IL+1.
224 *> If NZC = -1, then a workspace query is assumed; the
225 *> routine calculates the number of columns of the array Z that
226 *> are needed to hold the eigenvectors.
227 *> This value is returned as the first entry of the Z array, and
228 *> no error message related to NZC is issued by XERBLA.
229 *> \endverbatim
230 *>
231 *> \param[out] ISUPPZ
232 *> \verbatim
233 *> ISUPPZ is INTEGER ARRAY, dimension ( 2*max(1,M) )
234 *> The support of the eigenvectors in Z, i.e., the indices
235 *> indicating the nonzero elements in Z. The i-th computed eigenvector
236 *> is nonzero only in elements ISUPPZ( 2*i-1 ) through
237 *> ISUPPZ( 2*i ). This is relevant in the case when the matrix
238 *> is split. ISUPPZ is only accessed when JOBZ is 'V' and N > 0.
239 *> \endverbatim
240 *>
241 *> \param[in,out] TRYRAC
242 *> \verbatim
243 *> TRYRAC is LOGICAL
244 *> If TRYRAC.EQ..TRUE., indicates that the code should check whether
245 *> the tridiagonal matrix defines its eigenvalues to high relative
246 *> accuracy. If so, the code uses relative-accuracy preserving
247 *> algorithms that might be (a bit) slower depending on the matrix.
248 *> If the matrix does not define its eigenvalues to high relative
249 *> accuracy, the code can uses possibly faster algorithms.
250 *> If TRYRAC.EQ..FALSE., the code is not required to guarantee
251 *> relatively accurate eigenvalues and can use the fastest possible
252 *> techniques.
253 *> On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix
254 *> does not define its eigenvalues to high relative accuracy.
255 *> \endverbatim
256 *>
257 *> \param[out] WORK
258 *> \verbatim
259 *> WORK is DOUBLE PRECISION array, dimension (LWORK)
260 *> On exit, if INFO = 0, WORK(1) returns the optimal
261 *> (and minimal) LWORK.
262 *> \endverbatim
263 *>
264 *> \param[in] LWORK
265 *> \verbatim
266 *> LWORK is INTEGER
267 *> The dimension of the array WORK. LWORK >= max(1,18*N)
268 *> if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'.
269 *> If LWORK = -1, then a workspace query is assumed; the routine
270 *> only calculates the optimal size of the WORK array, returns
271 *> this value as the first entry of the WORK array, and no error
272 *> message related to LWORK is issued by XERBLA.
273 *> \endverbatim
274 *>
275 *> \param[out] IWORK
276 *> \verbatim
277 *> IWORK is INTEGER array, dimension (LIWORK)
278 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
279 *> \endverbatim
280 *>
281 *> \param[in] LIWORK
282 *> \verbatim
283 *> LIWORK is INTEGER
284 *> The dimension of the array IWORK. LIWORK >= max(1,10*N)
285 *> if the eigenvectors are desired, and LIWORK >= max(1,8*N)
286 *> if only the eigenvalues are to be computed.
287 *> If LIWORK = -1, then a workspace query is assumed; the
288 *> routine only calculates the optimal size of the IWORK array,
289 *> returns this value as the first entry of the IWORK array, and
290 *> no error message related to LIWORK is issued by XERBLA.
291 *> \endverbatim
292 *>
293 *> \param[out] INFO
294 *> \verbatim
295 *> INFO is INTEGER
296 *> On exit, INFO
297 *> = 0: successful exit
298 *> < 0: if INFO = -i, the i-th argument had an illegal value
299 *> > 0: if INFO = 1X, internal error in DLARRE,
300 *> if INFO = 2X, internal error in ZLARRV.
301 *> Here, the digit X = ABS( IINFO ) < 10, where IINFO is
302 *> the nonzero error code returned by DLARRE or
303 *> ZLARRV, respectively.
304 *> \endverbatim
305 *
306 * Authors:
307 * ========
308 *
309 *> \author Univ. of Tennessee
310 *> \author Univ. of California Berkeley
311 *> \author Univ. of Colorado Denver
312 *> \author NAG Ltd.
313 *
314 *> \date November 2013
315 *
316 *> \ingroup complex16OTHERcomputational
317 *
318 *> \par Contributors:
319 * ==================
320 *>
321 *> Beresford Parlett, University of California, Berkeley, USA \n
322 *> Jim Demmel, University of California, Berkeley, USA \n
323 *> Inderjit Dhillon, University of Texas, Austin, USA \n
324 *> Osni Marques, LBNL/NERSC, USA \n
325 *> Christof Voemel, University of California, Berkeley, USA \n
326 *
327 * =====================================================================
328  SUBROUTINE zstemr( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
329  $ m, w, z, ldz, nzc, isuppz, tryrac, work, lwork,
330  $ iwork, liwork, info )
331 *
332 * -- LAPACK computational routine (version 3.5.0) --
333 * -- LAPACK is a software package provided by Univ. of Tennessee, --
334 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
335 * November 2013
336 *
337 * .. Scalar Arguments ..
338  CHARACTER jobz, range
339  LOGICAL tryrac
340  INTEGER il, info, iu, ldz, nzc, liwork, lwork, m, n
341  DOUBLE PRECISION vl, vu
342 * ..
343 * .. Array Arguments ..
344  INTEGER isuppz( * ), iwork( * )
345  DOUBLE PRECISION d( * ), e( * ), w( * ), work( * )
346  COMPLEX*16 z( ldz, * )
347 * ..
348 *
349 * =====================================================================
350 *
351 * .. Parameters ..
352  DOUBLE PRECISION zero, one, four, minrgp
353  parameter( zero = 0.0d0, one = 1.0d0,
354  $ four = 4.0d0,
355  $ minrgp = 1.0d-3 )
356 * ..
357 * .. Local Scalars ..
358  LOGICAL alleig, indeig, lquery, valeig, wantz, zquery
359  INTEGER i, ibegin, iend, ifirst, iil, iindbl, iindw,
360  $ iindwk, iinfo, iinspl, iiu, ilast, in, indd,
361  $ inde2, inderr, indgp, indgrs, indwrk, itmp,
362  $ itmp2, j, jblk, jj, liwmin, lwmin, nsplit,
363  $ nzcmin, offset, wbegin, wend
364  DOUBLE PRECISION bignum, cs, eps, pivmin, r1, r2, rmax, rmin,
365  $ rtol1, rtol2, safmin, scale, smlnum, sn,
366  $ thresh, tmp, tnrm, wl, wu
367 * ..
368 * ..
369 * .. External Functions ..
370  LOGICAL lsame
371  DOUBLE PRECISION dlamch, dlanst
372  EXTERNAL lsame, dlamch, dlanst
373 * ..
374 * .. External Subroutines ..
375  EXTERNAL dcopy, dlae2, dlaev2, dlarrc, dlarre, dlarrj,
377 * ..
378 * .. Intrinsic Functions ..
379  INTRINSIC max, min, sqrt
380 
381 
382 * ..
383 * .. Executable Statements ..
384 *
385 * Test the input parameters.
386 *
387  wantz = lsame( jobz, 'V' )
388  alleig = lsame( range, 'A' )
389  valeig = lsame( range, 'V' )
390  indeig = lsame( range, 'I' )
391 *
392  lquery = ( ( lwork.EQ.-1 ).OR.( liwork.EQ.-1 ) )
393  zquery = ( nzc.EQ.-1 )
394 
395 * DSTEMR needs WORK of size 6*N, IWORK of size 3*N.
396 * In addition, DLARRE needs WORK of size 6*N, IWORK of size 5*N.
397 * Furthermore, ZLARRV needs WORK of size 12*N, IWORK of size 7*N.
398  IF( wantz ) THEN
399  lwmin = 18*n
400  liwmin = 10*n
401  ELSE
402 * need less workspace if only the eigenvalues are wanted
403  lwmin = 12*n
404  liwmin = 8*n
405  ENDIF
406 
407  wl = zero
408  wu = zero
409  iil = 0
410  iiu = 0
411  nsplit = 0
412 
413  IF( valeig ) THEN
414 * We do not reference VL, VU in the cases RANGE = 'I','A'
415 * The interval (WL, WU] contains all the wanted eigenvalues.
416 * It is either given by the user or computed in DLARRE.
417  wl = vl
418  wu = vu
419  ELSEIF( indeig ) THEN
420 * We do not reference IL, IU in the cases RANGE = 'V','A'
421  iil = il
422  iiu = iu
423  ENDIF
424 *
425  info = 0
426  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
427  info = -1
428  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
429  info = -2
430  ELSE IF( n.LT.0 ) THEN
431  info = -3
432  ELSE IF( valeig .AND. n.GT.0 .AND. wu.LE.wl ) THEN
433  info = -7
434  ELSE IF( indeig .AND. ( iil.LT.1 .OR. iil.GT.n ) ) THEN
435  info = -8
436  ELSE IF( indeig .AND. ( iiu.LT.iil .OR. iiu.GT.n ) ) THEN
437  info = -9
438  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
439  info = -13
440  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
441  info = -17
442  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
443  info = -19
444  END IF
445 *
446 * Get machine constants.
447 *
448  safmin = dlamch( 'Safe minimum' )
449  eps = dlamch( 'Precision' )
450  smlnum = safmin / eps
451  bignum = one / smlnum
452  rmin = sqrt( smlnum )
453  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
454 *
455  IF( info.EQ.0 ) THEN
456  work( 1 ) = lwmin
457  iwork( 1 ) = liwmin
458 *
459  IF( wantz .AND. alleig ) THEN
460  nzcmin = n
461  ELSE IF( wantz .AND. valeig ) THEN
462  CALL dlarrc( 'T', n, vl, vu, d, e, safmin,
463  $ nzcmin, itmp, itmp2, info )
464  ELSE IF( wantz .AND. indeig ) THEN
465  nzcmin = iiu-iil+1
466  ELSE
467 * WANTZ .EQ. FALSE.
468  nzcmin = 0
469  ENDIF
470  IF( zquery .AND. info.EQ.0 ) THEN
471  z( 1,1 ) = nzcmin
472  ELSE IF( nzc.LT.nzcmin .AND. .NOT.zquery ) THEN
473  info = -14
474  END IF
475  END IF
476 
477  IF( info.NE.0 ) THEN
478 *
479  CALL xerbla( 'ZSTEMR', -info )
480 *
481  RETURN
482  ELSE IF( lquery .OR. zquery ) THEN
483  RETURN
484  END IF
485 *
486 * Handle N = 0, 1, and 2 cases immediately
487 *
488  m = 0
489  IF( n.EQ.0 )
490  $ RETURN
491 *
492  IF( n.EQ.1 ) THEN
493  IF( alleig .OR. indeig ) THEN
494  m = 1
495  w( 1 ) = d( 1 )
496  ELSE
497  IF( wl.LT.d( 1 ) .AND. wu.GE.d( 1 ) ) THEN
498  m = 1
499  w( 1 ) = d( 1 )
500  END IF
501  END IF
502  IF( wantz.AND.(.NOT.zquery) ) THEN
503  z( 1, 1 ) = one
504  isuppz(1) = 1
505  isuppz(2) = 1
506  END IF
507  RETURN
508  END IF
509 *
510  IF( n.EQ.2 ) THEN
511  IF( .NOT.wantz ) THEN
512  CALL dlae2( d(1), e(1), d(2), r1, r2 )
513  ELSE IF( wantz.AND.(.NOT.zquery) ) THEN
514  CALL dlaev2( d(1), e(1), d(2), r1, r2, cs, sn )
515  END IF
516  IF( alleig.OR.
517  $ (valeig.AND.(r2.GT.wl).AND.
518  $ (r2.LE.wu)).OR.
519  $ (indeig.AND.(iil.EQ.1)) ) THEN
520  m = m+1
521  w( m ) = r2
522  IF( wantz.AND.(.NOT.zquery) ) THEN
523  z( 1, m ) = -sn
524  z( 2, m ) = cs
525 * Note: At most one of SN and CS can be zero.
526  IF (sn.NE.zero) THEN
527  IF (cs.NE.zero) THEN
528  isuppz(2*m-1) = 1
529  isuppz(2*m-1) = 2
530  ELSE
531  isuppz(2*m-1) = 1
532  isuppz(2*m-1) = 1
533  END IF
534  ELSE
535  isuppz(2*m-1) = 2
536  isuppz(2*m) = 2
537  END IF
538  ENDIF
539  ENDIF
540  IF( alleig.OR.
541  $ (valeig.AND.(r1.GT.wl).AND.
542  $ (r1.LE.wu)).OR.
543  $ (indeig.AND.(iiu.EQ.2)) ) THEN
544  m = m+1
545  w( m ) = r1
546  IF( wantz.AND.(.NOT.zquery) ) THEN
547  z( 1, m ) = cs
548  z( 2, m ) = sn
549 * Note: At most one of SN and CS can be zero.
550  IF (sn.NE.zero) THEN
551  IF (cs.NE.zero) THEN
552  isuppz(2*m-1) = 1
553  isuppz(2*m-1) = 2
554  ELSE
555  isuppz(2*m-1) = 1
556  isuppz(2*m-1) = 1
557  END IF
558  ELSE
559  isuppz(2*m-1) = 2
560  isuppz(2*m) = 2
561  END IF
562  ENDIF
563  ENDIF
564  ELSE
565 
566 * Continue with general N
567 
568  indgrs = 1
569  inderr = 2*n + 1
570  indgp = 3*n + 1
571  indd = 4*n + 1
572  inde2 = 5*n + 1
573  indwrk = 6*n + 1
574 *
575  iinspl = 1
576  iindbl = n + 1
577  iindw = 2*n + 1
578  iindwk = 3*n + 1
579 *
580 * Scale matrix to allowable range, if necessary.
581 * The allowable range is related to the PIVMIN parameter; see the
582 * comments in DLARRD. The preference for scaling small values
583 * up is heuristic; we expect users' matrices not to be close to the
584 * RMAX threshold.
585 *
586  scale = one
587  tnrm = dlanst( 'M', n, d, e )
588  IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
589  scale = rmin / tnrm
590  ELSE IF( tnrm.GT.rmax ) THEN
591  scale = rmax / tnrm
592  END IF
593  IF( scale.NE.one ) THEN
594  CALL dscal( n, scale, d, 1 )
595  CALL dscal( n-1, scale, e, 1 )
596  tnrm = tnrm*scale
597  IF( valeig ) THEN
598 * If eigenvalues in interval have to be found,
599 * scale (WL, WU] accordingly
600  wl = wl*scale
601  wu = wu*scale
602  ENDIF
603  END IF
604 *
605 * Compute the desired eigenvalues of the tridiagonal after splitting
606 * into smaller subblocks if the corresponding off-diagonal elements
607 * are small
608 * THRESH is the splitting parameter for DLARRE
609 * A negative THRESH forces the old splitting criterion based on the
610 * size of the off-diagonal. A positive THRESH switches to splitting
611 * which preserves relative accuracy.
612 *
613  IF( tryrac ) THEN
614 * Test whether the matrix warrants the more expensive relative approach.
615  CALL dlarrr( n, d, e, iinfo )
616  ELSE
617 * The user does not care about relative accurately eigenvalues
618  iinfo = -1
619  ENDIF
620 * Set the splitting criterion
621  IF (iinfo.EQ.0) THEN
622  thresh = eps
623  ELSE
624  thresh = -eps
625 * relative accuracy is desired but T does not guarantee it
626  tryrac = .false.
627  ENDIF
628 *
629  IF( tryrac ) THEN
630 * Copy original diagonal, needed to guarantee relative accuracy
631  CALL dcopy(n,d,1,work(indd),1)
632  ENDIF
633 * Store the squares of the offdiagonal values of T
634  DO 5 j = 1, n-1
635  work( inde2+j-1 ) = e(j)**2
636  5 CONTINUE
637 
638 * Set the tolerance parameters for bisection
639  IF( .NOT.wantz ) THEN
640 * DLARRE computes the eigenvalues to full precision.
641  rtol1 = four * eps
642  rtol2 = four * eps
643  ELSE
644 * DLARRE computes the eigenvalues to less than full precision.
645 * ZLARRV will refine the eigenvalue approximations, and we only
646 * need less accurate initial bisection in DLARRE.
647 * Note: these settings do only affect the subset case and DLARRE
648  rtol1 = sqrt(eps)
649  rtol2 = max( sqrt(eps)*5.0d-3, four * eps )
650  ENDIF
651  CALL dlarre( range, n, wl, wu, iil, iiu, d, e,
652  $ work(inde2), rtol1, rtol2, thresh, nsplit,
653  $ iwork( iinspl ), m, w, work( inderr ),
654  $ work( indgp ), iwork( iindbl ),
655  $ iwork( iindw ), work( indgrs ), pivmin,
656  $ work( indwrk ), iwork( iindwk ), iinfo )
657  IF( iinfo.NE.0 ) THEN
658  info = 10 + abs( iinfo )
659  RETURN
660  END IF
661 * Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired
662 * part of the spectrum. All desired eigenvalues are contained in
663 * (WL,WU]
664 
665 
666  IF( wantz ) THEN
667 *
668 * Compute the desired eigenvectors corresponding to the computed
669 * eigenvalues
670 *
671  CALL zlarrv( n, wl, wu, d, e,
672  $ pivmin, iwork( iinspl ), m,
673  $ 1, m, minrgp, rtol1, rtol2,
674  $ w, work( inderr ), work( indgp ), iwork( iindbl ),
675  $ iwork( iindw ), work( indgrs ), z, ldz,
676  $ isuppz, work( indwrk ), iwork( iindwk ), iinfo )
677  IF( iinfo.NE.0 ) THEN
678  info = 20 + abs( iinfo )
679  RETURN
680  END IF
681  ELSE
682 * DLARRE computes eigenvalues of the (shifted) root representation
683 * ZLARRV returns the eigenvalues of the unshifted matrix.
684 * However, if the eigenvectors are not desired by the user, we need
685 * to apply the corresponding shifts from DLARRE to obtain the
686 * eigenvalues of the original matrix.
687  DO 20 j = 1, m
688  itmp = iwork( iindbl+j-1 )
689  w( j ) = w( j ) + e( iwork( iinspl+itmp-1 ) )
690  20 CONTINUE
691  END IF
692 *
693 
694  IF ( tryrac ) THEN
695 * Refine computed eigenvalues so that they are relatively accurate
696 * with respect to the original matrix T.
697  ibegin = 1
698  wbegin = 1
699  DO 39 jblk = 1, iwork( iindbl+m-1 )
700  iend = iwork( iinspl+jblk-1 )
701  in = iend - ibegin + 1
702  wend = wbegin - 1
703 * check if any eigenvalues have to be refined in this block
704  36 CONTINUE
705  IF( wend.LT.m ) THEN
706  IF( iwork( iindbl+wend ).EQ.jblk ) THEN
707  wend = wend + 1
708  go to 36
709  END IF
710  END IF
711  IF( wend.LT.wbegin ) THEN
712  ibegin = iend + 1
713  go to 39
714  END IF
715 
716  offset = iwork(iindw+wbegin-1)-1
717  ifirst = iwork(iindw+wbegin-1)
718  ilast = iwork(iindw+wend-1)
719  rtol2 = four * eps
720  CALL dlarrj( in,
721  $ work(indd+ibegin-1), work(inde2+ibegin-1),
722  $ ifirst, ilast, rtol2, offset, w(wbegin),
723  $ work( inderr+wbegin-1 ),
724  $ work( indwrk ), iwork( iindwk ), pivmin,
725  $ tnrm, iinfo )
726  ibegin = iend + 1
727  wbegin = wend + 1
728  39 CONTINUE
729  ENDIF
730 *
731 * If matrix was scaled, then rescale eigenvalues appropriately.
732 *
733  IF( scale.NE.one ) THEN
734  CALL dscal( m, one / scale, w, 1 )
735  END IF
736  END IF
737 *
738 * If eigenvalues are not in increasing order, then sort them,
739 * possibly along with eigenvectors.
740 *
741  IF( nsplit.GT.1 .OR. n.EQ.2 ) THEN
742  IF( .NOT. wantz ) THEN
743  CALL dlasrt( 'I', m, w, iinfo )
744  IF( iinfo.NE.0 ) THEN
745  info = 3
746  RETURN
747  END IF
748  ELSE
749  DO 60 j = 1, m - 1
750  i = 0
751  tmp = w( j )
752  DO 50 jj = j + 1, m
753  IF( w( jj ).LT.tmp ) THEN
754  i = jj
755  tmp = w( jj )
756  END IF
757  50 CONTINUE
758  IF( i.NE.0 ) THEN
759  w( i ) = w( j )
760  w( j ) = tmp
761  IF( wantz ) THEN
762  CALL zswap( n, z( 1, i ), 1, z( 1, j ), 1 )
763  itmp = isuppz( 2*i-1 )
764  isuppz( 2*i-1 ) = isuppz( 2*j-1 )
765  isuppz( 2*j-1 ) = itmp
766  itmp = isuppz( 2*i )
767  isuppz( 2*i ) = isuppz( 2*j )
768  isuppz( 2*j ) = itmp
769  END IF
770  END IF
771  60 CONTINUE
772  END IF
773  ENDIF
774 *
775 *
776  work( 1 ) = lwmin
777  iwork( 1 ) = liwmin
778  RETURN
779 *
780 * End of ZSTEMR
781 *
782  END