LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
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 November 2013
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.5.0) --
316 * -- LAPACK is a software package provided by Univ. of Tennessee, --
317 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
318 * November 2013
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  nsplit = 0
395 
396  IF( valeig ) THEN
397 * We do not reference VL, VU in the cases RANGE = 'I','A'
398 * The interval (WL, WU] contains all the wanted eigenvalues.
399 * It is either given by the user or computed in DLARRE.
400  wl = vl
401  wu = vu
402  ELSEIF( indeig ) THEN
403 * We do not reference IL, IU in the cases RANGE = 'V','A'
404  iil = il
405  iiu = iu
406  ENDIF
407 *
408  info = 0
409  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
410  info = -1
411  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
412  info = -2
413  ELSE IF( n.LT.0 ) THEN
414  info = -3
415  ELSE IF( valeig .AND. n.GT.0 .AND. wu.LE.wl ) THEN
416  info = -7
417  ELSE IF( indeig .AND. ( iil.LT.1 .OR. iil.GT.n ) ) THEN
418  info = -8
419  ELSE IF( indeig .AND. ( iiu.LT.iil .OR. iiu.GT.n ) ) THEN
420  info = -9
421  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
422  info = -13
423  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
424  info = -17
425  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
426  info = -19
427  END IF
428 *
429 * Get machine constants.
430 *
431  safmin = dlamch( 'Safe minimum' )
432  eps = dlamch( 'Precision' )
433  smlnum = safmin / eps
434  bignum = one / smlnum
435  rmin = sqrt( smlnum )
436  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
437 *
438  IF( info.EQ.0 ) THEN
439  work( 1 ) = lwmin
440  iwork( 1 ) = liwmin
441 *
442  IF( wantz .AND. alleig ) THEN
443  nzcmin = n
444  ELSE IF( wantz .AND. valeig ) THEN
445  CALL dlarrc( 'T', n, vl, vu, d, e, safmin,
446  $ nzcmin, itmp, itmp2, info )
447  ELSE IF( wantz .AND. indeig ) THEN
448  nzcmin = iiu-iil+1
449  ELSE
450 * WANTZ .EQ. FALSE.
451  nzcmin = 0
452  ENDIF
453  IF( zquery .AND. info.EQ.0 ) THEN
454  z( 1,1 ) = nzcmin
455  ELSE IF( nzc.LT.nzcmin .AND. .NOT.zquery ) THEN
456  info = -14
457  END IF
458  END IF
459 
460  IF( info.NE.0 ) THEN
461 *
462  CALL xerbla( 'DSTEMR', -info )
463 *
464  RETURN
465  ELSE IF( lquery .OR. zquery ) THEN
466  RETURN
467  END IF
468 *
469 * Handle N = 0, 1, and 2 cases immediately
470 *
471  m = 0
472  IF( n.EQ.0 )
473  $ RETURN
474 *
475  IF( n.EQ.1 ) THEN
476  IF( alleig .OR. indeig ) THEN
477  m = 1
478  w( 1 ) = d( 1 )
479  ELSE
480  IF( wl.LT.d( 1 ) .AND. wu.GE.d( 1 ) ) THEN
481  m = 1
482  w( 1 ) = d( 1 )
483  END IF
484  END IF
485  IF( wantz.AND.(.NOT.zquery) ) THEN
486  z( 1, 1 ) = one
487  isuppz(1) = 1
488  isuppz(2) = 1
489  END IF
490  RETURN
491  END IF
492 *
493  IF( n.EQ.2 ) THEN
494  IF( .NOT.wantz ) THEN
495  CALL dlae2( d(1), e(1), d(2), r1, r2 )
496  ELSE IF( wantz.AND.(.NOT.zquery) ) THEN
497  CALL dlaev2( d(1), e(1), d(2), r1, r2, cs, sn )
498  END IF
499  IF( alleig.OR.
500  $ (valeig.AND.(r2.GT.wl).AND.
501  $ (r2.LE.wu)).OR.
502  $ (indeig.AND.(iil.EQ.1)) ) THEN
503  m = m+1
504  w( m ) = r2
505  IF( wantz.AND.(.NOT.zquery) ) THEN
506  z( 1, m ) = -sn
507  z( 2, m ) = cs
508 * Note: At most one of SN and CS can be zero.
509  IF (sn.NE.zero) THEN
510  IF (cs.NE.zero) THEN
511  isuppz(2*m-1) = 1
512  isuppz(2*m) = 2
513  ELSE
514  isuppz(2*m-1) = 1
515  isuppz(2*m) = 1
516  END IF
517  ELSE
518  isuppz(2*m-1) = 2
519  isuppz(2*m) = 2
520  END IF
521  ENDIF
522  ENDIF
523  IF( alleig.OR.
524  $ (valeig.AND.(r1.GT.wl).AND.
525  $ (r1.LE.wu)).OR.
526  $ (indeig.AND.(iiu.EQ.2)) ) THEN
527  m = m+1
528  w( m ) = r1
529  IF( wantz.AND.(.NOT.zquery) ) THEN
530  z( 1, m ) = cs
531  z( 2, m ) = sn
532 * Note: At most one of SN and CS can be zero.
533  IF (sn.NE.zero) THEN
534  IF (cs.NE.zero) THEN
535  isuppz(2*m-1) = 1
536  isuppz(2*m) = 2
537  ELSE
538  isuppz(2*m-1) = 1
539  isuppz(2*m) = 1
540  END IF
541  ELSE
542  isuppz(2*m-1) = 2
543  isuppz(2*m) = 2
544  END IF
545  ENDIF
546  ENDIF
547 
548  ELSE
549 
550 * Continue with general N
551 
552  indgrs = 1
553  inderr = 2*n + 1
554  indgp = 3*n + 1
555  indd = 4*n + 1
556  inde2 = 5*n + 1
557  indwrk = 6*n + 1
558 *
559  iinspl = 1
560  iindbl = n + 1
561  iindw = 2*n + 1
562  iindwk = 3*n + 1
563 *
564 * Scale matrix to allowable range, if necessary.
565 * The allowable range is related to the PIVMIN parameter; see the
566 * comments in DLARRD. The preference for scaling small values
567 * up is heuristic; we expect users' matrices not to be close to the
568 * RMAX threshold.
569 *
570  scale = one
571  tnrm = dlanst( 'M', n, d, e )
572  IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
573  scale = rmin / tnrm
574  ELSE IF( tnrm.GT.rmax ) THEN
575  scale = rmax / tnrm
576  END IF
577  IF( scale.NE.one ) THEN
578  CALL dscal( n, scale, d, 1 )
579  CALL dscal( n-1, scale, e, 1 )
580  tnrm = tnrm*scale
581  IF( valeig ) THEN
582 * If eigenvalues in interval have to be found,
583 * scale (WL, WU] accordingly
584  wl = wl*scale
585  wu = wu*scale
586  ENDIF
587  END IF
588 *
589 * Compute the desired eigenvalues of the tridiagonal after splitting
590 * into smaller subblocks if the corresponding off-diagonal elements
591 * are small
592 * THRESH is the splitting parameter for DLARRE
593 * A negative THRESH forces the old splitting criterion based on the
594 * size of the off-diagonal. A positive THRESH switches to splitting
595 * which preserves relative accuracy.
596 *
597  IF( tryrac ) THEN
598 * Test whether the matrix warrants the more expensive relative approach.
599  CALL dlarrr( n, d, e, iinfo )
600  ELSE
601 * The user does not care about relative accurately eigenvalues
602  iinfo = -1
603  ENDIF
604 * Set the splitting criterion
605  IF (iinfo.EQ.0) THEN
606  thresh = eps
607  ELSE
608  thresh = -eps
609 * relative accuracy is desired but T does not guarantee it
610  tryrac = .false.
611  ENDIF
612 *
613  IF( tryrac ) THEN
614 * Copy original diagonal, needed to guarantee relative accuracy
615  CALL dcopy(n,d,1,work(indd),1)
616  ENDIF
617 * Store the squares of the offdiagonal values of T
618  DO 5 j = 1, n-1
619  work( inde2+j-1 ) = e(j)**2
620  5 CONTINUE
621 
622 * Set the tolerance parameters for bisection
623  IF( .NOT.wantz ) THEN
624 * DLARRE computes the eigenvalues to full precision.
625  rtol1 = four * eps
626  rtol2 = four * eps
627  ELSE
628 * DLARRE computes the eigenvalues to less than full precision.
629 * DLARRV will refine the eigenvalue approximations, and we can
630 * need less accurate initial bisection in DLARRE.
631 * Note: these settings do only affect the subset case and DLARRE
632  rtol1 = sqrt(eps)
633  rtol2 = max( sqrt(eps)*5.0d-3, four * eps )
634  ENDIF
635  CALL dlarre( range, n, wl, wu, iil, iiu, d, e,
636  $ work(inde2), rtol1, rtol2, thresh, nsplit,
637  $ iwork( iinspl ), m, w, work( inderr ),
638  $ work( indgp ), iwork( iindbl ),
639  $ iwork( iindw ), work( indgrs ), pivmin,
640  $ work( indwrk ), iwork( iindwk ), iinfo )
641  IF( iinfo.NE.0 ) THEN
642  info = 10 + abs( iinfo )
643  RETURN
644  END IF
645 * Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired
646 * part of the spectrum. All desired eigenvalues are contained in
647 * (WL,WU]
648 
649 
650  IF( wantz ) THEN
651 *
652 * Compute the desired eigenvectors corresponding to the computed
653 * eigenvalues
654 *
655  CALL dlarrv( n, wl, wu, d, e,
656  $ pivmin, iwork( iinspl ), m,
657  $ 1, m, minrgp, rtol1, rtol2,
658  $ w, work( inderr ), work( indgp ), iwork( iindbl ),
659  $ iwork( iindw ), work( indgrs ), z, ldz,
660  $ isuppz, work( indwrk ), iwork( iindwk ), iinfo )
661  IF( iinfo.NE.0 ) THEN
662  info = 20 + abs( iinfo )
663  RETURN
664  END IF
665  ELSE
666 * DLARRE computes eigenvalues of the (shifted) root representation
667 * DLARRV returns the eigenvalues of the unshifted matrix.
668 * However, if the eigenvectors are not desired by the user, we need
669 * to apply the corresponding shifts from DLARRE to obtain the
670 * eigenvalues of the original matrix.
671  DO 20 j = 1, m
672  itmp = iwork( iindbl+j-1 )
673  w( j ) = w( j ) + e( iwork( iinspl+itmp-1 ) )
674  20 CONTINUE
675  END IF
676 *
677 
678  IF ( tryrac ) THEN
679 * Refine computed eigenvalues so that they are relatively accurate
680 * with respect to the original matrix T.
681  ibegin = 1
682  wbegin = 1
683  DO 39 jblk = 1, iwork( iindbl+m-1 )
684  iend = iwork( iinspl+jblk-1 )
685  in = iend - ibegin + 1
686  wend = wbegin - 1
687 * check if any eigenvalues have to be refined in this block
688  36 CONTINUE
689  IF( wend.LT.m ) THEN
690  IF( iwork( iindbl+wend ).EQ.jblk ) THEN
691  wend = wend + 1
692  go to 36
693  END IF
694  END IF
695  IF( wend.LT.wbegin ) THEN
696  ibegin = iend + 1
697  go to 39
698  END IF
699 
700  offset = iwork(iindw+wbegin-1)-1
701  ifirst = iwork(iindw+wbegin-1)
702  ilast = iwork(iindw+wend-1)
703  rtol2 = four * eps
704  CALL dlarrj( in,
705  $ work(indd+ibegin-1), work(inde2+ibegin-1),
706  $ ifirst, ilast, rtol2, offset, w(wbegin),
707  $ work( inderr+wbegin-1 ),
708  $ work( indwrk ), iwork( iindwk ), pivmin,
709  $ tnrm, iinfo )
710  ibegin = iend + 1
711  wbegin = wend + 1
712  39 CONTINUE
713  ENDIF
714 *
715 * If matrix was scaled, then rescale eigenvalues appropriately.
716 *
717  IF( scale.NE.one ) THEN
718  CALL dscal( m, one / scale, w, 1 )
719  END IF
720 
721  END IF
722 
723 *
724 * If eigenvalues are not in increasing order, then sort them,
725 * possibly along with eigenvectors.
726 *
727  IF( nsplit.GT.1 .OR. n.EQ.2 ) THEN
728  IF( .NOT. wantz ) THEN
729  CALL dlasrt( 'I', m, w, iinfo )
730  IF( iinfo.NE.0 ) THEN
731  info = 3
732  RETURN
733  END IF
734  ELSE
735  DO 60 j = 1, m - 1
736  i = 0
737  tmp = w( j )
738  DO 50 jj = j + 1, m
739  IF( w( jj ).LT.tmp ) THEN
740  i = jj
741  tmp = w( jj )
742  END IF
743  50 CONTINUE
744  IF( i.NE.0 ) THEN
745  w( i ) = w( j )
746  w( j ) = tmp
747  IF( wantz ) THEN
748  CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
749  itmp = isuppz( 2*i-1 )
750  isuppz( 2*i-1 ) = isuppz( 2*j-1 )
751  isuppz( 2*j-1 ) = itmp
752  itmp = isuppz( 2*i )
753  isuppz( 2*i ) = isuppz( 2*j )
754  isuppz( 2*j ) = itmp
755  END IF
756  END IF
757  60 CONTINUE
758  END IF
759  ENDIF
760 *
761 *
762  work( 1 ) = lwmin
763  iwork( 1 ) = liwmin
764  RETURN
765 *
766 * End of DSTEMR
767 *
768  END