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