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