LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 *>
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 REAL
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 REAL 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 REAL 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 = .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 = .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 REAL 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 SLARRE,
292 *> if INFO = 2X, internal error in SLARRV.
293 *> Here, the digit X = ABS( IINFO ) < 10, where IINFO is
294 *> the nonzero error code returned by SLARRE or
295 *> SLARRV, 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 *> \ingroup realOTHERcomputational
307 *
308 *> \par Contributors:
309 * ==================
310 *>
311 *> Beresford Parlett, University of California, Berkeley, USA \n
312 *> Jim Demmel, University of California, Berkeley, USA \n
313 *> Inderjit Dhillon, University of Texas, Austin, USA \n
314 *> Osni Marques, LBNL/NERSC, USA \n
315 *> Christof Voemel, University of California, Berkeley, USA
316 *
317 * =====================================================================
318  SUBROUTINE sstemr( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
319  $ M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK,
320  $ IWORK, LIWORK, INFO )
321 *
322 * -- LAPACK computational routine --
323 * -- LAPACK is a software package provided by Univ. of Tennessee, --
324 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
325 *
326 * .. Scalar Arguments ..
327  CHARACTER JOBZ, RANGE
328  LOGICAL TRYRAC
329  INTEGER IL, INFO, IU, LDZ, NZC, LIWORK, LWORK, M, N
330  REAL VL, VU
331 * ..
332 * .. Array Arguments ..
333  INTEGER ISUPPZ( * ), IWORK( * )
334  REAL D( * ), E( * ), W( * ), WORK( * )
335  REAL Z( LDZ, * )
336 * ..
337 *
338 * =====================================================================
339 *
340 * .. Parameters ..
341  REAL ZERO, ONE, FOUR, MINRGP
342  PARAMETER ( ZERO = 0.0e0, one = 1.0e0,
343  $ four = 4.0e0,
344  $ minrgp = 3.0e-3 )
345 * ..
346 * .. Local Scalars ..
347  LOGICAL ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY
348  INTEGER I, IBEGIN, IEND, IFIRST, IIL, IINDBL, IINDW,
349  $ iindwk, iinfo, iinspl, iiu, ilast, in, indd,
350  $ inde2, inderr, indgp, indgrs, indwrk, itmp,
351  $ itmp2, j, jblk, jj, liwmin, lwmin, nsplit,
352  $ nzcmin, offset, wbegin, wend
353  REAL BIGNUM, CS, EPS, PIVMIN, R1, R2, RMAX, RMIN,
354  $ RTOL1, RTOL2, SAFMIN, SCALE, SMLNUM, SN,
355  $ thresh, tmp, tnrm, wl, wu
356 * ..
357 * ..
358 * .. External Functions ..
359  LOGICAL LSAME
360  REAL SLAMCH, SLANST
361  EXTERNAL lsame, slamch, slanst
362 * ..
363 * .. External Subroutines ..
364  EXTERNAL scopy, slae2, slaev2, slarrc, slarre, slarrj,
366 * ..
367 * .. Intrinsic Functions ..
368  INTRINSIC max, min, sqrt
369 * ..
370 * .. Executable Statements ..
371 *
372 * Test the input parameters.
373 *
374  wantz = lsame( jobz, 'V' )
375  alleig = lsame( range, 'A' )
376  valeig = lsame( range, 'V' )
377  indeig = lsame( range, 'I' )
378 *
379  lquery = ( ( lwork.EQ.-1 ).OR.( liwork.EQ.-1 ) )
380  zquery = ( nzc.EQ.-1 )
381 
382 * SSTEMR needs WORK of size 6*N, IWORK of size 3*N.
383 * In addition, SLARRE needs WORK of size 6*N, IWORK of size 5*N.
384 * Furthermore, SLARRV needs WORK of size 12*N, IWORK of size 7*N.
385  IF( wantz ) THEN
386  lwmin = 18*n
387  liwmin = 10*n
388  ELSE
389 * need less workspace if only the eigenvalues are wanted
390  lwmin = 12*n
391  liwmin = 8*n
392  ENDIF
393 
394  wl = zero
395  wu = zero
396  iil = 0
397  iiu = 0
398  nsplit = 0
399 
400  IF( valeig ) THEN
401 * We do not reference VL, VU in the cases RANGE = 'I','A'
402 * The interval (WL, WU] contains all the wanted eigenvalues.
403 * It is either given by the user or computed in SLARRE.
404  wl = vl
405  wu = vu
406  ELSEIF( indeig ) THEN
407 * We do not reference IL, IU in the cases RANGE = 'V','A'
408  iil = il
409  iiu = iu
410  ENDIF
411 *
412  info = 0
413  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
414  info = -1
415  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
416  info = -2
417  ELSE IF( n.LT.0 ) THEN
418  info = -3
419  ELSE IF( valeig .AND. n.GT.0 .AND. wu.LE.wl ) THEN
420  info = -7
421  ELSE IF( indeig .AND. ( iil.LT.1 .OR. iil.GT.n ) ) THEN
422  info = -8
423  ELSE IF( indeig .AND. ( iiu.LT.iil .OR. iiu.GT.n ) ) THEN
424  info = -9
425  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
426  info = -13
427  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
428  info = -17
429  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
430  info = -19
431  END IF
432 *
433 * Get machine constants.
434 *
435  safmin = slamch( 'Safe minimum' )
436  eps = slamch( 'Precision' )
437  smlnum = safmin / eps
438  bignum = one / smlnum
439  rmin = sqrt( smlnum )
440  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
441 *
442  IF( info.EQ.0 ) THEN
443  work( 1 ) = lwmin
444  iwork( 1 ) = liwmin
445 *
446  IF( wantz .AND. alleig ) THEN
447  nzcmin = n
448  ELSE IF( wantz .AND. valeig ) THEN
449  CALL slarrc( 'T', n, vl, vu, d, e, safmin,
450  $ nzcmin, itmp, itmp2, info )
451  ELSE IF( wantz .AND. indeig ) THEN
452  nzcmin = iiu-iil+1
453  ELSE
454 * WANTZ .EQ. FALSE.
455  nzcmin = 0
456  ENDIF
457  IF( zquery .AND. info.EQ.0 ) THEN
458  z( 1,1 ) = nzcmin
459  ELSE IF( nzc.LT.nzcmin .AND. .NOT.zquery ) THEN
460  info = -14
461  END IF
462  END IF
463 
464  IF( info.NE.0 ) THEN
465 *
466  CALL xerbla( 'SSTEMR', -info )
467 *
468  RETURN
469  ELSE IF( lquery .OR. zquery ) THEN
470  RETURN
471  END IF
472 *
473 * Handle N = 0, 1, and 2 cases immediately
474 *
475  m = 0
476  IF( n.EQ.0 )
477  $ RETURN
478 *
479  IF( n.EQ.1 ) THEN
480  IF( alleig .OR. indeig ) THEN
481  m = 1
482  w( 1 ) = d( 1 )
483  ELSE
484  IF( wl.LT.d( 1 ) .AND. wu.GE.d( 1 ) ) THEN
485  m = 1
486  w( 1 ) = d( 1 )
487  END IF
488  END IF
489  IF( wantz.AND.(.NOT.zquery) ) THEN
490  z( 1, 1 ) = one
491  isuppz(1) = 1
492  isuppz(2) = 1
493  END IF
494  RETURN
495  END IF
496 *
497  IF( n.EQ.2 ) THEN
498  IF( .NOT.wantz ) THEN
499  CALL slae2( d(1), e(1), d(2), r1, r2 )
500  ELSE IF( wantz.AND.(.NOT.zquery) ) THEN
501  CALL slaev2( d(1), e(1), d(2), r1, r2, cs, sn )
502  END IF
503  IF( alleig.OR.
504  $ (valeig.AND.(r2.GT.wl).AND.
505  $ (r2.LE.wu)).OR.
506  $ (indeig.AND.(iil.EQ.1)) ) THEN
507  m = m+1
508  w( m ) = r2
509  IF( wantz.AND.(.NOT.zquery) ) THEN
510  z( 1, m ) = -sn
511  z( 2, m ) = cs
512 * Note: At most one of SN and CS can be zero.
513  IF (sn.NE.zero) THEN
514  IF (cs.NE.zero) THEN
515  isuppz(2*m-1) = 1
516  isuppz(2*m) = 2
517  ELSE
518  isuppz(2*m-1) = 1
519  isuppz(2*m) = 1
520  END IF
521  ELSE
522  isuppz(2*m-1) = 2
523  isuppz(2*m) = 2
524  END IF
525  ENDIF
526  ENDIF
527  IF( alleig.OR.
528  $ (valeig.AND.(r1.GT.wl).AND.
529  $ (r1.LE.wu)).OR.
530  $ (indeig.AND.(iiu.EQ.2)) ) THEN
531  m = m+1
532  w( m ) = r1
533  IF( wantz.AND.(.NOT.zquery) ) THEN
534  z( 1, m ) = cs
535  z( 2, m ) = sn
536 * Note: At most one of SN and CS can be zero.
537  IF (sn.NE.zero) THEN
538  IF (cs.NE.zero) THEN
539  isuppz(2*m-1) = 1
540  isuppz(2*m) = 2
541  ELSE
542  isuppz(2*m-1) = 1
543  isuppz(2*m) = 1
544  END IF
545  ELSE
546  isuppz(2*m-1) = 2
547  isuppz(2*m) = 2
548  END IF
549  ENDIF
550  ENDIF
551  ELSE
552 
553 * Continue with general N
554 
555  indgrs = 1
556  inderr = 2*n + 1
557  indgp = 3*n + 1
558  indd = 4*n + 1
559  inde2 = 5*n + 1
560  indwrk = 6*n + 1
561 *
562  iinspl = 1
563  iindbl = n + 1
564  iindw = 2*n + 1
565  iindwk = 3*n + 1
566 *
567 * Scale matrix to allowable range, if necessary.
568 * The allowable range is related to the PIVMIN parameter; see the
569 * comments in SLARRD. The preference for scaling small values
570 * up is heuristic; we expect users' matrices not to be close to the
571 * RMAX threshold.
572 *
573  scale = one
574  tnrm = slanst( 'M', n, d, e )
575  IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
576  scale = rmin / tnrm
577  ELSE IF( tnrm.GT.rmax ) THEN
578  scale = rmax / tnrm
579  END IF
580  IF( scale.NE.one ) THEN
581  CALL sscal( n, scale, d, 1 )
582  CALL sscal( n-1, scale, e, 1 )
583  tnrm = tnrm*scale
584  IF( valeig ) THEN
585 * If eigenvalues in interval have to be found,
586 * scale (WL, WU] accordingly
587  wl = wl*scale
588  wu = wu*scale
589  ENDIF
590  END IF
591 *
592 * Compute the desired eigenvalues of the tridiagonal after splitting
593 * into smaller subblocks if the corresponding off-diagonal elements
594 * are small
595 * THRESH is the splitting parameter for SLARRE
596 * A negative THRESH forces the old splitting criterion based on the
597 * size of the off-diagonal. A positive THRESH switches to splitting
598 * which preserves relative accuracy.
599 *
600  IF( tryrac ) THEN
601 * Test whether the matrix warrants the more expensive relative approach.
602  CALL slarrr( n, d, e, iinfo )
603  ELSE
604 * The user does not care about relative accurately eigenvalues
605  iinfo = -1
606  ENDIF
607 * Set the splitting criterion
608  IF (iinfo.EQ.0) THEN
609  thresh = eps
610  ELSE
611  thresh = -eps
612 * relative accuracy is desired but T does not guarantee it
613  tryrac = .false.
614  ENDIF
615 *
616  IF( tryrac ) THEN
617 * Copy original diagonal, needed to guarantee relative accuracy
618  CALL scopy(n,d,1,work(indd),1)
619  ENDIF
620 * Store the squares of the offdiagonal values of T
621  DO 5 j = 1, n-1
622  work( inde2+j-1 ) = e(j)**2
623  5 CONTINUE
624 
625 * Set the tolerance parameters for bisection
626  IF( .NOT.wantz ) THEN
627 * SLARRE computes the eigenvalues to full precision.
628  rtol1 = four * eps
629  rtol2 = four * eps
630  ELSE
631 * SLARRE computes the eigenvalues to less than full precision.
632 * SLARRV will refine the eigenvalue approximations, and we can
633 * need less accurate initial bisection in SLARRE.
634 * Note: these settings do only affect the subset case and SLARRE
635  rtol1 = max( sqrt(eps)*5.0e-2, four * eps )
636  rtol2 = max( sqrt(eps)*5.0e-3, four * eps )
637  ENDIF
638  CALL slarre( range, n, wl, wu, iil, iiu, d, e,
639  $ work(inde2), rtol1, rtol2, thresh, nsplit,
640  $ iwork( iinspl ), m, w, work( inderr ),
641  $ work( indgp ), iwork( iindbl ),
642  $ iwork( iindw ), work( indgrs ), pivmin,
643  $ work( indwrk ), iwork( iindwk ), iinfo )
644  IF( iinfo.NE.0 ) THEN
645  info = 10 + abs( iinfo )
646  RETURN
647  END IF
648 * Note that if RANGE .NE. 'V', SLARRE computes bounds on the desired
649 * part of the spectrum. All desired eigenvalues are contained in
650 * (WL,WU]
651 
652 
653  IF( wantz ) THEN
654 *
655 * Compute the desired eigenvectors corresponding to the computed
656 * eigenvalues
657 *
658  CALL slarrv( n, wl, wu, d, e,
659  $ pivmin, iwork( iinspl ), m,
660  $ 1, m, minrgp, rtol1, rtol2,
661  $ w, work( inderr ), work( indgp ), iwork( iindbl ),
662  $ iwork( iindw ), work( indgrs ), z, ldz,
663  $ isuppz, work( indwrk ), iwork( iindwk ), iinfo )
664  IF( iinfo.NE.0 ) THEN
665  info = 20 + abs( iinfo )
666  RETURN
667  END IF
668  ELSE
669 * SLARRE computes eigenvalues of the (shifted) root representation
670 * SLARRV returns the eigenvalues of the unshifted matrix.
671 * However, if the eigenvectors are not desired by the user, we need
672 * to apply the corresponding shifts from SLARRE to obtain the
673 * eigenvalues of the original matrix.
674  DO 20 j = 1, m
675  itmp = iwork( iindbl+j-1 )
676  w( j ) = w( j ) + e( iwork( iinspl+itmp-1 ) )
677  20 CONTINUE
678  END IF
679 *
680 
681  IF ( tryrac ) THEN
682 * Refine computed eigenvalues so that they are relatively accurate
683 * with respect to the original matrix T.
684  ibegin = 1
685  wbegin = 1
686  DO 39 jblk = 1, iwork( iindbl+m-1 )
687  iend = iwork( iinspl+jblk-1 )
688  in = iend - ibegin + 1
689  wend = wbegin - 1
690 * check if any eigenvalues have to be refined in this block
691  36 CONTINUE
692  IF( wend.LT.m ) THEN
693  IF( iwork( iindbl+wend ).EQ.jblk ) THEN
694  wend = wend + 1
695  GO TO 36
696  END IF
697  END IF
698  IF( wend.LT.wbegin ) THEN
699  ibegin = iend + 1
700  GO TO 39
701  END IF
702 
703  offset = iwork(iindw+wbegin-1)-1
704  ifirst = iwork(iindw+wbegin-1)
705  ilast = iwork(iindw+wend-1)
706  rtol2 = four * eps
707  CALL slarrj( in,
708  $ work(indd+ibegin-1), work(inde2+ibegin-1),
709  $ ifirst, ilast, rtol2, offset, w(wbegin),
710  $ work( inderr+wbegin-1 ),
711  $ work( indwrk ), iwork( iindwk ), pivmin,
712  $ tnrm, iinfo )
713  ibegin = iend + 1
714  wbegin = wend + 1
715  39 CONTINUE
716  ENDIF
717 *
718 * If matrix was scaled, then rescale eigenvalues appropriately.
719 *
720  IF( scale.NE.one ) THEN
721  CALL sscal( m, one / scale, w, 1 )
722  END IF
723  END IF
724 *
725 * If eigenvalues are not in increasing order, then sort them,
726 * possibly along with eigenvectors.
727 *
728  IF( nsplit.GT.1 .OR. n.EQ.2 ) THEN
729  IF( .NOT. wantz ) THEN
730  CALL slasrt( 'I', m, w, iinfo )
731  IF( iinfo.NE.0 ) THEN
732  info = 3
733  RETURN
734  END IF
735  ELSE
736  DO 60 j = 1, m - 1
737  i = 0
738  tmp = w( j )
739  DO 50 jj = j + 1, m
740  IF( w( jj ).LT.tmp ) THEN
741  i = jj
742  tmp = w( jj )
743  END IF
744  50 CONTINUE
745  IF( i.NE.0 ) THEN
746  w( i ) = w( j )
747  w( j ) = tmp
748  IF( wantz ) THEN
749  CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
750  itmp = isuppz( 2*i-1 )
751  isuppz( 2*i-1 ) = isuppz( 2*j-1 )
752  isuppz( 2*j-1 ) = itmp
753  itmp = isuppz( 2*i )
754  isuppz( 2*i ) = isuppz( 2*j )
755  isuppz( 2*j ) = itmp
756  END IF
757  END IF
758  60 CONTINUE
759  END IF
760  ENDIF
761 *
762 *
763  work( 1 ) = lwmin
764  iwork( 1 ) = liwmin
765  RETURN
766 *
767 * End of SSTEMR
768 *
769  END
subroutine slarrr(N, D, E, INFO)
SLARRR performs tests to decide whether the symmetric tridiagonal matrix T warrants expensive computa...
Definition: slarrr.f:94
subroutine slarrc(JOBT, N, VL, VU, D, E, PIVMIN, EIGCNT, LCNT, RCNT, INFO)
SLARRC computes the number of eigenvalues of the symmetric tridiagonal matrix.
Definition: slarrc.f:137
subroutine slarre(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)
SLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduce...
Definition: slarre.f:305
subroutine slarrj(N, D, E2, IFIRST, ILAST, RTOL, OFFSET, W, WERR, WORK, IWORK, PIVMIN, SPDIAM, INFO)
SLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T.
Definition: slarrj.f:168
subroutine slae2(A, B, C, RT1, RT2)
SLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition: slae2.f:102
subroutine slaev2(A, B, C, RT1, RT2, CS1, SN1)
SLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.
Definition: slaev2.f:120
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine slasrt(ID, N, D, INFO)
SLASRT sorts numbers in increasing or decreasing order.
Definition: slasrt.f:88
subroutine slarrv(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)
SLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues ...
Definition: slarrv.f:292
subroutine sstemr(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, M, W, Z, LDZ, NZC, ISUPPZ, TRYRAC, WORK, LWORK, IWORK, LIWORK, INFO)
SSTEMR
Definition: sstemr.f:321
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79