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