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