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