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