LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zstemr.f
Go to the documentation of this file.
1 *> \brief \b ZSTEMR
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZSTEMR + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zstemr.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zstemr.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zstemr.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZSTEMR( 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 * DOUBLE PRECISION VL, VU
30 * ..
31 * .. Array Arguments ..
32 * INTEGER ISUPPZ( * ), IWORK( * )
33 * DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
34 * COMPLEX*16 Z( LDZ, * )
35 * ..
36 *
37 *
38 *> \par Purpose:
39 * =============
40 *>
41 *> \verbatim
42 *>
43 *> ZSTEMR 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.ZSTEMR 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 *> ZSTEMR accepts complex workspace to facilitate interoperability
109 *> with ZUNMTR or ZUPMTR.
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
156 *> \endverbatim
157 *>
158 *> \param[in] VU
159 *> \verbatim
160 *> VU is DOUBLE PRECISION
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 DOUBLE PRECISION 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*16 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 DOUBLE PRECISION 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 DLARRE,
300 *> if INFO = 2X, internal error in ZLARRV.
301 *> Here, the digit X = ABS( IINFO ) < 10, where IINFO is
302 *> the nonzero error code returned by DLARRE or
303 *> ZLARRV, 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 September 2012
315 *
316 *> \ingroup complex16OTHERcomputational
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 \n
326 *
327 * =====================================================================
328  SUBROUTINE zstemr( 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.4.2) --
333 * -- LAPACK is a software package provided by Univ. of Tennessee, --
334 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
335 * September 2012
336 *
337 * .. Scalar Arguments ..
338  CHARACTER jobz, range
339  LOGICAL tryrac
340  INTEGER il, info, iu, ldz, nzc, liwork, lwork, m, n
341  DOUBLE PRECISION vl, vu
342 * ..
343 * .. Array Arguments ..
344  INTEGER isuppz( * ), iwork( * )
345  DOUBLE PRECISION d( * ), e( * ), w( * ), work( * )
346  COMPLEX*16 z( ldz, * )
347 * ..
348 *
349 * =====================================================================
350 *
351 * .. Parameters ..
352  DOUBLE PRECISION zero, one, four, minrgp
353  parameter( zero = 0.0d0, one = 1.0d0,
354  $ four = 4.0d0,
355  $ minrgp = 1.0d-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  DOUBLE PRECISION 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  DOUBLE PRECISION dlamch, dlanst
372  EXTERNAL lsame, dlamch, dlanst
373 * ..
374 * .. External Subroutines ..
375  EXTERNAL dcopy, dlae2, dlaev2, dlarrc, dlarre, dlarrj,
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 * DSTEMR needs WORK of size 6*N, IWORK of size 3*N.
396 * In addition, DLARRE needs WORK of size 6*N, IWORK of size 5*N.
397 * Furthermore, ZLARRV 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 
412  IF( valeig ) THEN
413 * We do not reference VL, VU in the cases RANGE = 'I','A'
414 * The interval (WL, WU] contains all the wanted eigenvalues.
415 * It is either given by the user or computed in DLARRE.
416  wl = vl
417  wu = vu
418  elseif( indeig ) THEN
419 * We do not reference IL, IU in the cases RANGE = 'V','A'
420  iil = il
421  iiu = iu
422  ENDIF
423 *
424  info = 0
425  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
426  info = -1
427  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
428  info = -2
429  ELSE IF( n.LT.0 ) THEN
430  info = -3
431  ELSE IF( valeig .AND. n.GT.0 .AND. wu.LE.wl ) THEN
432  info = -7
433  ELSE IF( indeig .AND. ( iil.LT.1 .OR. iil.GT.n ) ) THEN
434  info = -8
435  ELSE IF( indeig .AND. ( iiu.LT.iil .OR. iiu.GT.n ) ) THEN
436  info = -9
437  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
438  info = -13
439  ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
440  info = -17
441  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
442  info = -19
443  END IF
444 *
445 * Get machine constants.
446 *
447  safmin = dlamch( 'Safe minimum' )
448  eps = dlamch( 'Precision' )
449  smlnum = safmin / eps
450  bignum = one / smlnum
451  rmin = sqrt( smlnum )
452  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
453 *
454  IF( info.EQ.0 ) THEN
455  work( 1 ) = lwmin
456  iwork( 1 ) = liwmin
457 *
458  IF( wantz .AND. alleig ) THEN
459  nzcmin = n
460  ELSE IF( wantz .AND. valeig ) THEN
461  CALL dlarrc( 'T', n, vl, vu, d, e, safmin,
462  $ nzcmin, itmp, itmp2, info )
463  ELSE IF( wantz .AND. indeig ) THEN
464  nzcmin = iiu-iil+1
465  ELSE
466 * WANTZ .EQ. FALSE.
467  nzcmin = 0
468  ENDIF
469  IF( zquery .AND. info.EQ.0 ) THEN
470  z( 1,1 ) = nzcmin
471  ELSE IF( nzc.LT.nzcmin .AND. .NOT.zquery ) THEN
472  info = -14
473  END IF
474  END IF
475 
476  IF( info.NE.0 ) THEN
477 *
478  CALL xerbla( 'ZSTEMR', -info )
479 *
480  return
481  ELSE IF( lquery .OR. zquery ) THEN
482  return
483  END IF
484 *
485 * Handle N = 0, 1, and 2 cases immediately
486 *
487  m = 0
488  IF( n.EQ.0 )
489  $ return
490 *
491  IF( n.EQ.1 ) THEN
492  IF( alleig .OR. indeig ) THEN
493  m = 1
494  w( 1 ) = d( 1 )
495  ELSE
496  IF( wl.LT.d( 1 ) .AND. wu.GE.d( 1 ) ) THEN
497  m = 1
498  w( 1 ) = d( 1 )
499  END IF
500  END IF
501  IF( wantz.AND.(.NOT.zquery) ) THEN
502  z( 1, 1 ) = one
503  isuppz(1) = 1
504  isuppz(2) = 1
505  END IF
506  return
507  END IF
508 *
509  IF( n.EQ.2 ) THEN
510  IF( .NOT.wantz ) THEN
511  CALL dlae2( d(1), e(1), d(2), r1, r2 )
512  ELSE IF( wantz.AND.(.NOT.zquery) ) THEN
513  CALL dlaev2( d(1), e(1), d(2), r1, r2, cs, sn )
514  END IF
515  IF( alleig.OR.
516  $ (valeig.AND.(r2.GT.wl).AND.
517  $ (r2.LE.wu)).OR.
518  $ (indeig.AND.(iil.EQ.1)) ) THEN
519  m = m+1
520  w( m ) = r2
521  IF( wantz.AND.(.NOT.zquery) ) THEN
522  z( 1, m ) = -sn
523  z( 2, m ) = cs
524 * Note: At most one of SN and CS can be zero.
525  IF (sn.NE.zero) THEN
526  IF (cs.NE.zero) THEN
527  isuppz(2*m-1) = 1
528  isuppz(2*m-1) = 2
529  ELSE
530  isuppz(2*m-1) = 1
531  isuppz(2*m-1) = 1
532  END IF
533  ELSE
534  isuppz(2*m-1) = 2
535  isuppz(2*m) = 2
536  END IF
537  ENDIF
538  ENDIF
539  IF( alleig.OR.
540  $ (valeig.AND.(r1.GT.wl).AND.
541  $ (r1.LE.wu)).OR.
542  $ (indeig.AND.(iiu.EQ.2)) ) THEN
543  m = m+1
544  w( m ) = r1
545  IF( wantz.AND.(.NOT.zquery) ) THEN
546  z( 1, m ) = cs
547  z( 2, m ) = sn
548 * Note: At most one of SN and CS can be zero.
549  IF (sn.NE.zero) THEN
550  IF (cs.NE.zero) THEN
551  isuppz(2*m-1) = 1
552  isuppz(2*m-1) = 2
553  ELSE
554  isuppz(2*m-1) = 1
555  isuppz(2*m-1) = 1
556  END IF
557  ELSE
558  isuppz(2*m-1) = 2
559  isuppz(2*m) = 2
560  END IF
561  ENDIF
562  ENDIF
563  ELSE
564 
565 * Continue with general N
566 
567  indgrs = 1
568  inderr = 2*n + 1
569  indgp = 3*n + 1
570  indd = 4*n + 1
571  inde2 = 5*n + 1
572  indwrk = 6*n + 1
573 *
574  iinspl = 1
575  iindbl = n + 1
576  iindw = 2*n + 1
577  iindwk = 3*n + 1
578 *
579 * Scale matrix to allowable range, if necessary.
580 * The allowable range is related to the PIVMIN parameter; see the
581 * comments in DLARRD. The preference for scaling small values
582 * up is heuristic; we expect users' matrices not to be close to the
583 * RMAX threshold.
584 *
585  scale = one
586  tnrm = dlanst( 'M', n, d, e )
587  IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
588  scale = rmin / tnrm
589  ELSE IF( tnrm.GT.rmax ) THEN
590  scale = rmax / tnrm
591  END IF
592  IF( scale.NE.one ) THEN
593  CALL dscal( n, scale, d, 1 )
594  CALL dscal( n-1, scale, e, 1 )
595  tnrm = tnrm*scale
596  IF( valeig ) THEN
597 * If eigenvalues in interval have to be found,
598 * scale (WL, WU] accordingly
599  wl = wl*scale
600  wu = wu*scale
601  ENDIF
602  END IF
603 *
604 * Compute the desired eigenvalues of the tridiagonal after splitting
605 * into smaller subblocks if the corresponding off-diagonal elements
606 * are small
607 * THRESH is the splitting parameter for DLARRE
608 * A negative THRESH forces the old splitting criterion based on the
609 * size of the off-diagonal. A positive THRESH switches to splitting
610 * which preserves relative accuracy.
611 *
612  IF( tryrac ) THEN
613 * Test whether the matrix warrants the more expensive relative approach.
614  CALL dlarrr( n, d, e, iinfo )
615  ELSE
616 * The user does not care about relative accurately eigenvalues
617  iinfo = -1
618  ENDIF
619 * Set the splitting criterion
620  IF (iinfo.EQ.0) THEN
621  thresh = eps
622  ELSE
623  thresh = -eps
624 * relative accuracy is desired but T does not guarantee it
625  tryrac = .false.
626  ENDIF
627 *
628  IF( tryrac ) THEN
629 * Copy original diagonal, needed to guarantee relative accuracy
630  CALL dcopy(n,d,1,work(indd),1)
631  ENDIF
632 * Store the squares of the offdiagonal values of T
633  DO 5 j = 1, n-1
634  work( inde2+j-1 ) = e(j)**2
635  5 continue
636 
637 * Set the tolerance parameters for bisection
638  IF( .NOT.wantz ) THEN
639 * DLARRE computes the eigenvalues to full precision.
640  rtol1 = four * eps
641  rtol2 = four * eps
642  ELSE
643 * DLARRE computes the eigenvalues to less than full precision.
644 * ZLARRV will refine the eigenvalue approximations, and we only
645 * need less accurate initial bisection in DLARRE.
646 * Note: these settings do only affect the subset case and DLARRE
647  rtol1 = sqrt(eps)
648  rtol2 = max( sqrt(eps)*5.0d-3, four * eps )
649  ENDIF
650  CALL dlarre( range, n, wl, wu, iil, iiu, d, e,
651  $ work(inde2), rtol1, rtol2, thresh, nsplit,
652  $ iwork( iinspl ), m, w, work( inderr ),
653  $ work( indgp ), iwork( iindbl ),
654  $ iwork( iindw ), work( indgrs ), pivmin,
655  $ work( indwrk ), iwork( iindwk ), iinfo )
656  IF( iinfo.NE.0 ) THEN
657  info = 10 + abs( iinfo )
658  return
659  END IF
660 * Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired
661 * part of the spectrum. All desired eigenvalues are contained in
662 * (WL,WU]
663 
664 
665  IF( wantz ) THEN
666 *
667 * Compute the desired eigenvectors corresponding to the computed
668 * eigenvalues
669 *
670  CALL zlarrv( n, wl, wu, d, e,
671  $ pivmin, iwork( iinspl ), m,
672  $ 1, m, minrgp, rtol1, rtol2,
673  $ w, work( inderr ), work( indgp ), iwork( iindbl ),
674  $ iwork( iindw ), work( indgrs ), z, ldz,
675  $ isuppz, work( indwrk ), iwork( iindwk ), iinfo )
676  IF( iinfo.NE.0 ) THEN
677  info = 20 + abs( iinfo )
678  return
679  END IF
680  ELSE
681 * DLARRE computes eigenvalues of the (shifted) root representation
682 * ZLARRV returns the eigenvalues of the unshifted matrix.
683 * However, if the eigenvectors are not desired by the user, we need
684 * to apply the corresponding shifts from DLARRE to obtain the
685 * eigenvalues of the original matrix.
686  DO 20 j = 1, m
687  itmp = iwork( iindbl+j-1 )
688  w( j ) = w( j ) + e( iwork( iinspl+itmp-1 ) )
689  20 continue
690  END IF
691 *
692 
693  IF ( tryrac ) THEN
694 * Refine computed eigenvalues so that they are relatively accurate
695 * with respect to the original matrix T.
696  ibegin = 1
697  wbegin = 1
698  DO 39 jblk = 1, iwork( iindbl+m-1 )
699  iend = iwork( iinspl+jblk-1 )
700  in = iend - ibegin + 1
701  wend = wbegin - 1
702 * check if any eigenvalues have to be refined in this block
703  36 continue
704  IF( wend.LT.m ) THEN
705  IF( iwork( iindbl+wend ).EQ.jblk ) THEN
706  wend = wend + 1
707  go to 36
708  END IF
709  END IF
710  IF( wend.LT.wbegin ) THEN
711  ibegin = iend + 1
712  go to 39
713  END IF
714 
715  offset = iwork(iindw+wbegin-1)-1
716  ifirst = iwork(iindw+wbegin-1)
717  ilast = iwork(iindw+wend-1)
718  rtol2 = four * eps
719  CALL dlarrj( in,
720  $ work(indd+ibegin-1), work(inde2+ibegin-1),
721  $ ifirst, ilast, rtol2, offset, w(wbegin),
722  $ work( inderr+wbegin-1 ),
723  $ work( indwrk ), iwork( iindwk ), pivmin,
724  $ tnrm, iinfo )
725  ibegin = iend + 1
726  wbegin = wend + 1
727  39 continue
728  ENDIF
729 *
730 * If matrix was scaled, then rescale eigenvalues appropriately.
731 *
732  IF( scale.NE.one ) THEN
733  CALL dscal( m, one / scale, w, 1 )
734  END IF
735  END IF
736 *
737 * If eigenvalues are not in increasing order, then sort them,
738 * possibly along with eigenvectors.
739 *
740  IF( nsplit.GT.1 .OR. n.EQ.2 ) THEN
741  IF( .NOT. wantz ) THEN
742  CALL dlasrt( 'I', m, w, iinfo )
743  IF( iinfo.NE.0 ) THEN
744  info = 3
745  return
746  END IF
747  ELSE
748  DO 60 j = 1, m - 1
749  i = 0
750  tmp = w( j )
751  DO 50 jj = j + 1, m
752  IF( w( jj ).LT.tmp ) THEN
753  i = jj
754  tmp = w( jj )
755  END IF
756  50 continue
757  IF( i.NE.0 ) THEN
758  w( i ) = w( j )
759  w( j ) = tmp
760  IF( wantz ) THEN
761  CALL zswap( n, z( 1, i ), 1, z( 1, j ), 1 )
762  itmp = isuppz( 2*i-1 )
763  isuppz( 2*i-1 ) = isuppz( 2*j-1 )
764  isuppz( 2*j-1 ) = itmp
765  itmp = isuppz( 2*i )
766  isuppz( 2*i ) = isuppz( 2*j )
767  isuppz( 2*j ) = itmp
768  END IF
769  END IF
770  60 continue
771  END IF
772  ENDIF
773 *
774 *
775  work( 1 ) = lwmin
776  iwork( 1 ) = liwmin
777  return
778 *
779 * End of ZSTEMR
780 *
781  END