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