LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
sstemr.f
Go to the documentation of this file.
1 *> \brief \b SSTEMR
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SSTEMR + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sstemr.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sstemr.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sstemr.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SSTEMR( 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 * REAL VL, VU
30 * ..
31 * .. Array Arguments ..
32 * INTEGER ISUPPZ( * ), IWORK( * )
33 * REAL D( * ), E( * ), W( * ), WORK( * )
34 * REAL Z( LDZ, * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> SSTEMR 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.SSTEMR 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 REAL 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 REAL 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 REAL
139 *> \endverbatim
140 *>
141 *> \param[in] VU
142 *> \verbatim
143 *> VU is REAL
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 REAL 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 REAL 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 REAL 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 SLARRE,
283 *> if INFO = 2X, internal error in SLARRV.
284 *> Here, the digit X = ABS( IINFO ) < 10, where IINFO is
285 *> the nonzero error code returned by SLARRE or
286 *> SLARRV, 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 realOTHERcomputational
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 sstemr( 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  REAL vl, vu
325 * ..
326 * .. Array Arguments ..
327  INTEGER isuppz( * ), iwork( * )
328  REAL d( * ), e( * ), w( * ), work( * )
329  REAL z( ldz, * )
330 * ..
331 *
332 * =====================================================================
333 *
334 * .. Parameters ..
335  REAL zero, one, four, minrgp
336  parameter( zero = 0.0e0, one = 1.0e0,
337  $ four = 4.0e0,
338  $ minrgp = 3.0e-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  REAL 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  REAL slamch, slanst
355  EXTERNAL lsame, slamch, slanst
356 * ..
357 * .. External Subroutines ..
358  EXTERNAL scopy, slae2, slaev2, slarrc, slarre, slarrj,
360 * ..
361 * .. Intrinsic Functions ..
362  INTRINSIC max, min, sqrt
363 * ..
364 * .. Executable Statements ..
365 *
366 * Test the input parameters.
367 *
368  wantz = lsame( jobz, 'V' )
369  alleig = lsame( range, 'A' )
370  valeig = lsame( range, 'V' )
371  indeig = lsame( range, 'I' )
372 *
373  lquery = ( ( lwork.EQ.-1 ).OR.( liwork.EQ.-1 ) )
374  zquery = ( nzc.EQ.-1 )
375 
376 * SSTEMR needs WORK of size 6*N, IWORK of size 3*N.
377 * In addition, SLARRE needs WORK of size 6*N, IWORK of size 5*N.
378 * Furthermore, SLARRV needs WORK of size 12*N, IWORK of size 7*N.
379  IF( wantz ) THEN
380  lwmin = 18*n
381  liwmin = 10*n
382  ELSE
383 * need less workspace if only the eigenvalues are wanted
384  lwmin = 12*n
385  liwmin = 8*n
386  ENDIF
387 
388  wl = zero
389  wu = zero
390  iil = 0
391  iiu = 0
392 
393  IF( valeig ) THEN
394 * We do not reference VL, VU in the cases RANGE = 'I','A'
395 * The interval (WL, WU] contains all the wanted eigenvalues.
396 * It is either given by the user or computed in SLARRE.
397  wl = vl
398  wu = vu
399  elseif( indeig ) THEN
400 * We do not reference IL, IU in the cases RANGE = 'V','A'
401  iil = il
402  iiu = iu
403  ENDIF
404 *
405  info = 0
406  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
407  info = -1
408  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
409  info = -2
410  ELSE IF( n.LT.0 ) THEN
411  info = -3
412  ELSE IF( valeig .AND. n.GT.0 .AND. wu.LE.wl ) THEN
413  info = -7
414  ELSE IF( indeig .AND. ( iil.LT.1 .OR. iil.GT.n ) ) THEN
415  info = -8
416  ELSE IF( indeig .AND. ( iiu.LT.iil .OR. iiu.GT.n ) ) THEN
417  info = -9
418  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
419  info = -13
420  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
421  info = -17
422  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
423  info = -19
424  END IF
425 *
426 * Get machine constants.
427 *
428  safmin = slamch( 'Safe minimum' )
429  eps = slamch( 'Precision' )
430  smlnum = safmin / eps
431  bignum = one / smlnum
432  rmin = sqrt( smlnum )
433  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
434 *
435  IF( info.EQ.0 ) THEN
436  work( 1 ) = lwmin
437  iwork( 1 ) = liwmin
438 *
439  IF( wantz .AND. alleig ) THEN
440  nzcmin = n
441  ELSE IF( wantz .AND. valeig ) THEN
442  CALL slarrc( 'T', n, vl, vu, d, e, safmin,
443  $ nzcmin, itmp, itmp2, info )
444  ELSE IF( wantz .AND. indeig ) THEN
445  nzcmin = iiu-iil+1
446  ELSE
447 * WANTZ .EQ. FALSE.
448  nzcmin = 0
449  ENDIF
450  IF( zquery .AND. info.EQ.0 ) THEN
451  z( 1,1 ) = nzcmin
452  ELSE IF( nzc.LT.nzcmin .AND. .NOT.zquery ) THEN
453  info = -14
454  END IF
455  END IF
456 
457  IF( info.NE.0 ) THEN
458 *
459  CALL xerbla( 'SSTEMR', -info )
460 *
461  return
462  ELSE IF( lquery .OR. zquery ) THEN
463  return
464  END IF
465 *
466 * Handle N = 0, 1, and 2 cases immediately
467 *
468  m = 0
469  IF( n.EQ.0 )
470  $ return
471 *
472  IF( n.EQ.1 ) THEN
473  IF( alleig .OR. indeig ) THEN
474  m = 1
475  w( 1 ) = d( 1 )
476  ELSE
477  IF( wl.LT.d( 1 ) .AND. wu.GE.d( 1 ) ) THEN
478  m = 1
479  w( 1 ) = d( 1 )
480  END IF
481  END IF
482  IF( wantz.AND.(.NOT.zquery) ) THEN
483  z( 1, 1 ) = one
484  isuppz(1) = 1
485  isuppz(2) = 1
486  END IF
487  return
488  END IF
489 *
490  IF( n.EQ.2 ) THEN
491  IF( .NOT.wantz ) THEN
492  CALL slae2( d(1), e(1), d(2), r1, r2 )
493  ELSE IF( wantz.AND.(.NOT.zquery) ) THEN
494  CALL slaev2( d(1), e(1), d(2), r1, r2, cs, sn )
495  END IF
496  IF( alleig.OR.
497  $ (valeig.AND.(r2.GT.wl).AND.
498  $ (r2.LE.wu)).OR.
499  $ (indeig.AND.(iil.EQ.1)) ) THEN
500  m = m+1
501  w( m ) = r2
502  IF( wantz.AND.(.NOT.zquery) ) THEN
503  z( 1, m ) = -sn
504  z( 2, m ) = cs
505 * Note: At most one of SN and CS can be zero.
506  IF (sn.NE.zero) THEN
507  IF (cs.NE.zero) THEN
508  isuppz(2*m-1) = 1
509  isuppz(2*m) = 2
510  ELSE
511  isuppz(2*m-1) = 1
512  isuppz(2*m) = 1
513  END IF
514  ELSE
515  isuppz(2*m-1) = 2
516  isuppz(2*m) = 2
517  END IF
518  ENDIF
519  ENDIF
520  IF( alleig.OR.
521  $ (valeig.AND.(r1.GT.wl).AND.
522  $ (r1.LE.wu)).OR.
523  $ (indeig.AND.(iiu.EQ.2)) ) THEN
524  m = m+1
525  w( m ) = r1
526  IF( wantz.AND.(.NOT.zquery) ) THEN
527  z( 1, m ) = cs
528  z( 2, m ) = sn
529 * Note: At most one of SN and CS can be zero.
530  IF (sn.NE.zero) THEN
531  IF (cs.NE.zero) THEN
532  isuppz(2*m-1) = 1
533  isuppz(2*m) = 2
534  ELSE
535  isuppz(2*m-1) = 1
536  isuppz(2*m) = 1
537  END IF
538  ELSE
539  isuppz(2*m-1) = 2
540  isuppz(2*m) = 2
541  END IF
542  ENDIF
543  ENDIF
544  ELSE
545 
546 * Continue with general N
547 
548  indgrs = 1
549  inderr = 2*n + 1
550  indgp = 3*n + 1
551  indd = 4*n + 1
552  inde2 = 5*n + 1
553  indwrk = 6*n + 1
554 *
555  iinspl = 1
556  iindbl = n + 1
557  iindw = 2*n + 1
558  iindwk = 3*n + 1
559 *
560 * Scale matrix to allowable range, if necessary.
561 * The allowable range is related to the PIVMIN parameter; see the
562 * comments in SLARRD. The preference for scaling small values
563 * up is heuristic; we expect users' matrices not to be close to the
564 * RMAX threshold.
565 *
566  scale = one
567  tnrm = slanst( 'M', n, d, e )
568  IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
569  scale = rmin / tnrm
570  ELSE IF( tnrm.GT.rmax ) THEN
571  scale = rmax / tnrm
572  END IF
573  IF( scale.NE.one ) THEN
574  CALL sscal( n, scale, d, 1 )
575  CALL sscal( n-1, scale, e, 1 )
576  tnrm = tnrm*scale
577  IF( valeig ) THEN
578 * If eigenvalues in interval have to be found,
579 * scale (WL, WU] accordingly
580  wl = wl*scale
581  wu = wu*scale
582  ENDIF
583  END IF
584 *
585 * Compute the desired eigenvalues of the tridiagonal after splitting
586 * into smaller subblocks if the corresponding off-diagonal elements
587 * are small
588 * THRESH is the splitting parameter for SLARRE
589 * A negative THRESH forces the old splitting criterion based on the
590 * size of the off-diagonal. A positive THRESH switches to splitting
591 * which preserves relative accuracy.
592 *
593  IF( tryrac ) THEN
594 * Test whether the matrix warrants the more expensive relative approach.
595  CALL slarrr( n, d, e, iinfo )
596  ELSE
597 * The user does not care about relative accurately eigenvalues
598  iinfo = -1
599  ENDIF
600 * Set the splitting criterion
601  IF (iinfo.EQ.0) THEN
602  thresh = eps
603  ELSE
604  thresh = -eps
605 * relative accuracy is desired but T does not guarantee it
606  tryrac = .false.
607  ENDIF
608 *
609  IF( tryrac ) THEN
610 * Copy original diagonal, needed to guarantee relative accuracy
611  CALL scopy(n,d,1,work(indd),1)
612  ENDIF
613 * Store the squares of the offdiagonal values of T
614  DO 5 j = 1, n-1
615  work( inde2+j-1 ) = e(j)**2
616  5 continue
617 
618 * Set the tolerance parameters for bisection
619  IF( .NOT.wantz ) THEN
620 * SLARRE computes the eigenvalues to full precision.
621  rtol1 = four * eps
622  rtol2 = four * eps
623  ELSE
624 * SLARRE computes the eigenvalues to less than full precision.
625 * SLARRV will refine the eigenvalue approximations, and we can
626 * need less accurate initial bisection in SLARRE.
627 * Note: these settings do only affect the subset case and SLARRE
628  rtol1 = max( sqrt(eps)*5.0e-2, four * eps )
629  rtol2 = max( sqrt(eps)*5.0e-3, four * eps )
630  ENDIF
631  CALL slarre( range, n, wl, wu, iil, iiu, d, e,
632  $ work(inde2), rtol1, rtol2, thresh, nsplit,
633  $ iwork( iinspl ), m, w, work( inderr ),
634  $ work( indgp ), iwork( iindbl ),
635  $ iwork( iindw ), work( indgrs ), pivmin,
636  $ work( indwrk ), iwork( iindwk ), iinfo )
637  IF( iinfo.NE.0 ) THEN
638  info = 10 + abs( iinfo )
639  return
640  END IF
641 * Note that if RANGE .NE. 'V', SLARRE computes bounds on the desired
642 * part of the spectrum. All desired eigenvalues are contained in
643 * (WL,WU]
644 
645 
646  IF( wantz ) THEN
647 *
648 * Compute the desired eigenvectors corresponding to the computed
649 * eigenvalues
650 *
651  CALL slarrv( n, wl, wu, d, e,
652  $ pivmin, iwork( iinspl ), m,
653  $ 1, m, minrgp, rtol1, rtol2,
654  $ w, work( inderr ), work( indgp ), iwork( iindbl ),
655  $ iwork( iindw ), work( indgrs ), z, ldz,
656  $ isuppz, work( indwrk ), iwork( iindwk ), iinfo )
657  IF( iinfo.NE.0 ) THEN
658  info = 20 + abs( iinfo )
659  return
660  END IF
661  ELSE
662 * SLARRE computes eigenvalues of the (shifted) root representation
663 * SLARRV returns the eigenvalues of the unshifted matrix.
664 * However, if the eigenvectors are not desired by the user, we need
665 * to apply the corresponding shifts from SLARRE to obtain the
666 * eigenvalues of the original matrix.
667  DO 20 j = 1, m
668  itmp = iwork( iindbl+j-1 )
669  w( j ) = w( j ) + e( iwork( iinspl+itmp-1 ) )
670  20 continue
671  END IF
672 *
673 
674  IF ( tryrac ) THEN
675 * Refine computed eigenvalues so that they are relatively accurate
676 * with respect to the original matrix T.
677  ibegin = 1
678  wbegin = 1
679  DO 39 jblk = 1, iwork( iindbl+m-1 )
680  iend = iwork( iinspl+jblk-1 )
681  in = iend - ibegin + 1
682  wend = wbegin - 1
683 * check if any eigenvalues have to be refined in this block
684  36 continue
685  IF( wend.LT.m ) THEN
686  IF( iwork( iindbl+wend ).EQ.jblk ) THEN
687  wend = wend + 1
688  go to 36
689  END IF
690  END IF
691  IF( wend.LT.wbegin ) THEN
692  ibegin = iend + 1
693  go to 39
694  END IF
695 
696  offset = iwork(iindw+wbegin-1)-1
697  ifirst = iwork(iindw+wbegin-1)
698  ilast = iwork(iindw+wend-1)
699  rtol2 = four * eps
700  CALL slarrj( in,
701  $ work(indd+ibegin-1), work(inde2+ibegin-1),
702  $ ifirst, ilast, rtol2, offset, w(wbegin),
703  $ work( inderr+wbegin-1 ),
704  $ work( indwrk ), iwork( iindwk ), pivmin,
705  $ tnrm, iinfo )
706  ibegin = iend + 1
707  wbegin = wend + 1
708  39 continue
709  ENDIF
710 *
711 * If matrix was scaled, then rescale eigenvalues appropriately.
712 *
713  IF( scale.NE.one ) THEN
714  CALL sscal( m, one / scale, w, 1 )
715  END IF
716  END IF
717 *
718 * If eigenvalues are not in increasing order, then sort them,
719 * possibly along with eigenvectors.
720 *
721  IF( nsplit.GT.1 .OR. n.EQ.2 ) THEN
722  IF( .NOT. wantz ) THEN
723  CALL slasrt( 'I', m, w, iinfo )
724  IF( iinfo.NE.0 ) THEN
725  info = 3
726  return
727  END IF
728  ELSE
729  DO 60 j = 1, m - 1
730  i = 0
731  tmp = w( j )
732  DO 50 jj = j + 1, m
733  IF( w( jj ).LT.tmp ) THEN
734  i = jj
735  tmp = w( jj )
736  END IF
737  50 continue
738  IF( i.NE.0 ) THEN
739  w( i ) = w( j )
740  w( j ) = tmp
741  IF( wantz ) THEN
742  CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
743  itmp = isuppz( 2*i-1 )
744  isuppz( 2*i-1 ) = isuppz( 2*j-1 )
745  isuppz( 2*j-1 ) = itmp
746  itmp = isuppz( 2*i )
747  isuppz( 2*i ) = isuppz( 2*j )
748  isuppz( 2*j ) = itmp
749  END IF
750  END IF
751  60 continue
752  END IF
753  ENDIF
754 *
755 *
756  work( 1 ) = lwmin
757  iwork( 1 ) = liwmin
758  return
759 *
760 * End of SSTEMR
761 *
762  END