LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
slarrd.f
Go to the documentation of this file.
1 *> \brief \b SLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SLARRD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarrd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarrd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarrd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLARRD( RANGE, ORDER, N, VL, VU, IL, IU, GERS,
22 * RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
23 * M, W, WERR, WL, WU, IBLOCK, INDEXW,
24 * WORK, IWORK, INFO )
25 *
26 * .. Scalar Arguments ..
27 * CHARACTER ORDER, RANGE
28 * INTEGER IL, INFO, IU, M, N, NSPLIT
29 * REAL PIVMIN, RELTOL, VL, VU, WL, WU
30 * ..
31 * .. Array Arguments ..
32 * INTEGER IBLOCK( * ), INDEXW( * ),
33 * $ ISPLIT( * ), IWORK( * )
34 * REAL D( * ), E( * ), E2( * ),
35 * $ GERS( * ), W( * ), WERR( * ), WORK( * )
36 * ..
37 *
38 *
39 *> \par Purpose:
40 * =============
41 *>
42 *> \verbatim
43 *>
44 *> SLARRD computes the eigenvalues of a symmetric tridiagonal
45 *> matrix T to suitable accuracy. This is an auxiliary code to be
46 *> called from SSTEMR.
47 *> The user may ask for all eigenvalues, all eigenvalues
48 *> in the half-open interval (VL, VU], or the IL-th through IU-th
49 *> eigenvalues.
50 *>
51 *> To avoid overflow, the matrix must be scaled so that its
52 *> largest element is no greater than overflow**(1/2) * underflow**(1/4) in absolute value, and for greatest
53 *> accuracy, it should not be much smaller than that.
54 *>
55 *> See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal
56 *> Matrix", Report CS41, Computer Science Dept., Stanford
57 *> University, July 21, 1966.
58 *> \endverbatim
59 *
60 * Arguments:
61 * ==========
62 *
63 *> \param[in] RANGE
64 *> \verbatim
65 *> RANGE is CHARACTER*1
66 *> = 'A': ("All") all eigenvalues will be found.
67 *> = 'V': ("Value") all eigenvalues in the half-open interval
68 *> (VL, VU] will be found.
69 *> = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
70 *> entire matrix) will be found.
71 *> \endverbatim
72 *>
73 *> \param[in] ORDER
74 *> \verbatim
75 *> ORDER is CHARACTER*1
76 *> = 'B': ("By Block") the eigenvalues will be grouped by
77 *> split-off block (see IBLOCK, ISPLIT) and
78 *> ordered from smallest to largest within
79 *> the block.
80 *> = 'E': ("Entire matrix")
81 *> the eigenvalues for the entire matrix
82 *> will be ordered from smallest to
83 *> largest.
84 *> \endverbatim
85 *>
86 *> \param[in] N
87 *> \verbatim
88 *> N is INTEGER
89 *> The order of the tridiagonal matrix T. N >= 0.
90 *> \endverbatim
91 *>
92 *> \param[in] VL
93 *> \verbatim
94 *> VL is REAL
95 *> If RANGE='V', the lower bound of the interval to
96 *> be searched for eigenvalues. Eigenvalues less than or equal
97 *> to VL, or greater than VU, will not be returned. VL < VU.
98 *> Not referenced if RANGE = 'A' or 'I'.
99 *> \endverbatim
100 *>
101 *> \param[in] VU
102 *> \verbatim
103 *> VU is REAL
104 *> If RANGE='V', the upper bound of the interval to
105 *> be searched for eigenvalues. Eigenvalues less than or equal
106 *> to VL, or greater than VU, will not be returned. VL < VU.
107 *> Not referenced if RANGE = 'A' or 'I'.
108 *> \endverbatim
109 *>
110 *> \param[in] IL
111 *> \verbatim
112 *> IL is INTEGER
113 *> If RANGE='I', the index of the
114 *> smallest eigenvalue to be returned.
115 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
116 *> Not referenced if RANGE = 'A' or 'V'.
117 *> \endverbatim
118 *>
119 *> \param[in] IU
120 *> \verbatim
121 *> IU is INTEGER
122 *> If RANGE='I', the index of the
123 *> largest eigenvalue to be returned.
124 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
125 *> Not referenced if RANGE = 'A' or 'V'.
126 *> \endverbatim
127 *>
128 *> \param[in] GERS
129 *> \verbatim
130 *> GERS is REAL array, dimension (2*N)
131 *> The N Gerschgorin intervals (the i-th Gerschgorin interval
132 *> is (GERS(2*i-1), GERS(2*i)).
133 *> \endverbatim
134 *>
135 *> \param[in] RELTOL
136 *> \verbatim
137 *> RELTOL is REAL
138 *> The minimum relative width of an interval. When an interval
139 *> is narrower than RELTOL times the larger (in
140 *> magnitude) endpoint, then it is considered to be
141 *> sufficiently small, i.e., converged. Note: this should
142 *> always be at least radix*machine epsilon.
143 *> \endverbatim
144 *>
145 *> \param[in] D
146 *> \verbatim
147 *> D is REAL array, dimension (N)
148 *> The n diagonal elements of the tridiagonal matrix T.
149 *> \endverbatim
150 *>
151 *> \param[in] E
152 *> \verbatim
153 *> E is REAL array, dimension (N-1)
154 *> The (n-1) off-diagonal elements of the tridiagonal matrix T.
155 *> \endverbatim
156 *>
157 *> \param[in] E2
158 *> \verbatim
159 *> E2 is REAL array, dimension (N-1)
160 *> The (n-1) squared off-diagonal elements of the tridiagonal matrix T.
161 *> \endverbatim
162 *>
163 *> \param[in] PIVMIN
164 *> \verbatim
165 *> PIVMIN is REAL
166 *> The minimum pivot allowed in the Sturm sequence for T.
167 *> \endverbatim
168 *>
169 *> \param[in] NSPLIT
170 *> \verbatim
171 *> NSPLIT is INTEGER
172 *> The number of diagonal blocks in the matrix T.
173 *> 1 <= NSPLIT <= N.
174 *> \endverbatim
175 *>
176 *> \param[in] ISPLIT
177 *> \verbatim
178 *> ISPLIT is INTEGER array, dimension (N)
179 *> The splitting points, at which T breaks up into submatrices.
180 *> The first submatrix consists of rows/columns 1 to ISPLIT(1),
181 *> the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
182 *> etc., and the NSPLIT-th consists of rows/columns
183 *> ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
184 *> (Only the first NSPLIT elements will actually be used, but
185 *> since the user cannot know a priori what value NSPLIT will
186 *> have, N words must be reserved for ISPLIT.)
187 *> \endverbatim
188 *>
189 *> \param[out] M
190 *> \verbatim
191 *> M is INTEGER
192 *> The actual number of eigenvalues found. 0 <= M <= N.
193 *> (See also the description of INFO=2,3.)
194 *> \endverbatim
195 *>
196 *> \param[out] W
197 *> \verbatim
198 *> W is REAL array, dimension (N)
199 *> On exit, the first M elements of W will contain the
200 *> eigenvalue approximations. SLARRD computes an interval
201 *> I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue
202 *> approximation is given as the interval midpoint
203 *> W(j)= ( a_j + b_j)/2. The corresponding error is bounded by
204 *> WERR(j) = abs( a_j - b_j)/2
205 *> \endverbatim
206 *>
207 *> \param[out] WERR
208 *> \verbatim
209 *> WERR is REAL array, dimension (N)
210 *> The error bound on the corresponding eigenvalue approximation
211 *> in W.
212 *> \endverbatim
213 *>
214 *> \param[out] WL
215 *> \verbatim
216 *> WL is REAL
217 *> \endverbatim
218 *>
219 *> \param[out] WU
220 *> \verbatim
221 *> WU is REAL
222 *> The interval (WL, WU] contains all the wanted eigenvalues.
223 *> If RANGE='V', then WL=VL and WU=VU.
224 *> If RANGE='A', then WL and WU are the global Gerschgorin bounds
225 *> on the spectrum.
226 *> If RANGE='I', then WL and WU are computed by SLAEBZ from the
227 *> index range specified.
228 *> \endverbatim
229 *>
230 *> \param[out] IBLOCK
231 *> \verbatim
232 *> IBLOCK is INTEGER array, dimension (N)
233 *> At each row/column j where E(j) is zero or small, the
234 *> matrix T is considered to split into a block diagonal
235 *> matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which
236 *> block (from 1 to the number of blocks) the eigenvalue W(i)
237 *> belongs. (SLARRD may use the remaining N-M elements as
238 *> workspace.)
239 *> \endverbatim
240 *>
241 *> \param[out] INDEXW
242 *> \verbatim
243 *> INDEXW is INTEGER array, dimension (N)
244 *> The indices of the eigenvalues within each block (submatrix);
245 *> for example, INDEXW(i)= j and IBLOCK(i)=k imply that the
246 *> i-th eigenvalue W(i) is the j-th eigenvalue in block k.
247 *> \endverbatim
248 *>
249 *> \param[out] WORK
250 *> \verbatim
251 *> WORK is REAL array, dimension (4*N)
252 *> \endverbatim
253 *>
254 *> \param[out] IWORK
255 *> \verbatim
256 *> IWORK is INTEGER array, dimension (3*N)
257 *> \endverbatim
258 *>
259 *> \param[out] INFO
260 *> \verbatim
261 *> INFO is INTEGER
262 *> = 0: successful exit
263 *> < 0: if INFO = -i, the i-th argument had an illegal value
264 *> > 0: some or all of the eigenvalues failed to converge or
265 *> were not computed:
266 *> =1 or 3: Bisection failed to converge for some
267 *> eigenvalues; these eigenvalues are flagged by a
268 *> negative block number. The effect is that the
269 *> eigenvalues may not be as accurate as the
270 *> absolute and relative tolerances. This is
271 *> generally caused by unexpectedly inaccurate
272 *> arithmetic.
273 *> =2 or 3: RANGE='I' only: Not all of the eigenvalues
274 *> IL:IU were found.
275 *> Effect: M < IU+1-IL
276 *> Cause: non-monotonic arithmetic, causing the
277 *> Sturm sequence to be non-monotonic.
278 *> Cure: recalculate, using RANGE='A', and pick
279 *> out eigenvalues IL:IU. In some cases,
280 *> increasing the PARAMETER "FUDGE" may
281 *> make things work.
282 *> = 4: RANGE='I', and the Gershgorin interval
283 *> initially used was too small. No eigenvalues
284 *> were computed.
285 *> Probable cause: your machine has sloppy
286 *> floating-point arithmetic.
287 *> Cure: Increase the PARAMETER "FUDGE",
288 *> recompile, and try again.
289 *> \endverbatim
290 *
291 *> \par Internal Parameters:
292 * =========================
293 *>
294 *> \verbatim
295 *> FUDGE REAL, default = 2
296 *> A "fudge factor" to widen the Gershgorin intervals. Ideally,
297 *> a value of 1 should work, but on machines with sloppy
298 *> arithmetic, this needs to be larger. The default for
299 *> publicly released versions should be large enough to handle
300 *> the worst machine around. Note that this has no effect
301 *> on accuracy of the solution.
302 *> \endverbatim
303 *>
304 *> \par Contributors:
305 * ==================
306 *>
307 *> W. Kahan, University of California, Berkeley, USA \n
308 *> Beresford Parlett, University of California, Berkeley, USA \n
309 *> Jim Demmel, University of California, Berkeley, USA \n
310 *> Inderjit Dhillon, University of Texas, Austin, USA \n
311 *> Osni Marques, LBNL/NERSC, USA \n
312 *> Christof Voemel, University of California, Berkeley, USA \n
313 *
314 * Authors:
315 * ========
316 *
317 *> \author Univ. of Tennessee
318 *> \author Univ. of California Berkeley
319 *> \author Univ. of Colorado Denver
320 *> \author NAG Ltd.
321 *
322 *> \date June 2016
323 *
324 *> \ingroup OTHERauxiliary
325 *
326 * =====================================================================
327  SUBROUTINE slarrd( RANGE, ORDER, N, VL, VU, IL, IU, GERS,
328  $ RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT,
329  $ M, W, WERR, WL, WU, IBLOCK, INDEXW,
330  $ WORK, IWORK, INFO )
331 *
332 * -- LAPACK auxiliary routine (version 3.7.1) --
333 * -- LAPACK is a software package provided by Univ. of Tennessee, --
334 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
335 * June 2016
336 *
337 * .. Scalar Arguments ..
338  CHARACTER ORDER, RANGE
339  INTEGER IL, INFO, IU, M, N, NSPLIT
340  REAL PIVMIN, RELTOL, VL, VU, WL, WU
341 * ..
342 * .. Array Arguments ..
343  INTEGER IBLOCK( * ), INDEXW( * ),
344  $ ISPLIT( * ), IWORK( * )
345  REAL D( * ), E( * ), E2( * ),
346  $ gers( * ), w( * ), werr( * ), work( * )
347 * ..
348 *
349 * =====================================================================
350 *
351 * .. Parameters ..
352  REAL ZERO, ONE, TWO, HALF, FUDGE
353  PARAMETER ( ZERO = 0.0e0, one = 1.0e0,
354  $ two = 2.0e0, half = one/two,
355  $ fudge = two )
356  INTEGER ALLRNG, VALRNG, INDRNG
357  PARAMETER ( ALLRNG = 1, valrng = 2, indrng = 3 )
358 * ..
359 * .. Local Scalars ..
360  LOGICAL NCNVRG, TOOFEW
361  INTEGER I, IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO,
362  $ IM, IN, IOFF, IOUT, IRANGE, ITMAX, ITMP1,
363  $ itmp2, iw, iwoff, j, jblk, jdisc, je, jee, nb,
364  $ nwl, nwu
365  REAL ATOLI, EPS, GL, GU, RTOLI, TMP1, TMP2,
366  $ TNORM, UFLOW, WKILL, WLU, WUL
367 
368 * ..
369 * .. Local Arrays ..
370  INTEGER IDUMMA( 1 )
371 * ..
372 * .. External Functions ..
373  LOGICAL LSAME
374  INTEGER ILAENV
375  REAL SLAMCH
376  EXTERNAL lsame, ilaenv, slamch
377 * ..
378 * .. External Subroutines ..
379  EXTERNAL slaebz
380 * ..
381 * .. Intrinsic Functions ..
382  INTRINSIC abs, int, log, max, min
383 * ..
384 * .. Executable Statements ..
385 *
386  info = 0
387 *
388 * Quick return if possible
389 *
390  IF( n.LE.0 ) THEN
391  RETURN
392  END IF
393 *
394 * Decode RANGE
395 *
396  IF( lsame( range, 'A' ) ) THEN
397  irange = allrng
398  ELSE IF( lsame( range, 'V' ) ) THEN
399  irange = valrng
400  ELSE IF( lsame( range, 'I' ) ) THEN
401  irange = indrng
402  ELSE
403  irange = 0
404  END IF
405 *
406 * Check for Errors
407 *
408  IF( irange.LE.0 ) THEN
409  info = -1
410  ELSE IF( .NOT.(lsame(order,'B').OR.lsame(order,'E')) ) THEN
411  info = -2
412  ELSE IF( n.LT.0 ) THEN
413  info = -3
414  ELSE IF( irange.EQ.valrng ) THEN
415  IF( vl.GE.vu )
416  $ info = -5
417  ELSE IF( irange.EQ.indrng .AND.
418  $ ( il.LT.1 .OR. il.GT.max( 1, n ) ) ) THEN
419  info = -6
420  ELSE IF( irange.EQ.indrng .AND.
421  $ ( iu.LT.min( n, il ) .OR. iu.GT.n ) ) THEN
422  info = -7
423  END IF
424 *
425  IF( info.NE.0 ) THEN
426  RETURN
427  END IF
428 
429 * Initialize error flags
430  info = 0
431  ncnvrg = .false.
432  toofew = .false.
433 
434 * Quick return if possible
435  m = 0
436  IF( n.EQ.0 ) RETURN
437 
438 * Simplification:
439  IF( irange.EQ.indrng .AND. il.EQ.1 .AND. iu.EQ.n ) irange = 1
440 
441 * Get machine constants
442  eps = slamch( 'P' )
443  uflow = slamch( 'U' )
444 
445 
446 * Special Case when N=1
447 * Treat case of 1x1 matrix for quick return
448  IF( n.EQ.1 ) THEN
449  IF( (irange.EQ.allrng).OR.
450  $ ((irange.EQ.valrng).AND.(d(1).GT.vl).AND.(d(1).LE.vu)).OR.
451  $ ((irange.EQ.indrng).AND.(il.EQ.1).AND.(iu.EQ.1)) ) THEN
452  m = 1
453  w(1) = d(1)
454 * The computation error of the eigenvalue is zero
455  werr(1) = zero
456  iblock( 1 ) = 1
457  indexw( 1 ) = 1
458  ENDIF
459  RETURN
460  END IF
461 
462 * NB is the minimum vector length for vector bisection, or 0
463 * if only scalar is to be done.
464  nb = ilaenv( 1, 'SSTEBZ', ' ', n, -1, -1, -1 )
465  IF( nb.LE.1 ) nb = 0
466 
467 * Find global spectral radius
468  gl = d(1)
469  gu = d(1)
470  DO 5 i = 1,n
471  gl = min( gl, gers( 2*i - 1))
472  gu = max( gu, gers(2*i) )
473  5 CONTINUE
474 * Compute global Gerschgorin bounds and spectral diameter
475  tnorm = max( abs( gl ), abs( gu ) )
476  gl = gl - fudge*tnorm*eps*n - fudge*two*pivmin
477  gu = gu + fudge*tnorm*eps*n + fudge*two*pivmin
478 * [JAN/28/2009] remove the line below since SPDIAM variable not use
479 * SPDIAM = GU - GL
480 * Input arguments for SLAEBZ:
481 * The relative tolerance. An interval (a,b] lies within
482 * "relative tolerance" if b-a < RELTOL*max(|a|,|b|),
483  rtoli = reltol
484 * Set the absolute tolerance for interval convergence to zero to force
485 * interval convergence based on relative size of the interval.
486 * This is dangerous because intervals might not converge when RELTOL is
487 * small. But at least a very small number should be selected so that for
488 * strongly graded matrices, the code can get relatively accurate
489 * eigenvalues.
490  atoli = fudge*two*uflow + fudge*two*pivmin
491 
492  IF( irange.EQ.indrng ) THEN
493 
494 * RANGE='I': Compute an interval containing eigenvalues
495 * IL through IU. The initial interval [GL,GU] from the global
496 * Gerschgorin bounds GL and GU is refined by SLAEBZ.
497  itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /
498  $ log( two ) ) + 2
499  work( n+1 ) = gl
500  work( n+2 ) = gl
501  work( n+3 ) = gu
502  work( n+4 ) = gu
503  work( n+5 ) = gl
504  work( n+6 ) = gu
505  iwork( 1 ) = -1
506  iwork( 2 ) = -1
507  iwork( 3 ) = n + 1
508  iwork( 4 ) = n + 1
509  iwork( 5 ) = il - 1
510  iwork( 6 ) = iu
511 *
512  CALL slaebz( 3, itmax, n, 2, 2, nb, atoli, rtoli, pivmin,
513  $ d, e, e2, iwork( 5 ), work( n+1 ), work( n+5 ), iout,
514  $ iwork, w, iblock, iinfo )
515  IF( iinfo .NE. 0 ) THEN
516  info = iinfo
517  RETURN
518  END IF
519 * On exit, output intervals may not be ordered by ascending negcount
520  IF( iwork( 6 ).EQ.iu ) THEN
521  wl = work( n+1 )
522  wlu = work( n+3 )
523  nwl = iwork( 1 )
524  wu = work( n+4 )
525  wul = work( n+2 )
526  nwu = iwork( 4 )
527  ELSE
528  wl = work( n+2 )
529  wlu = work( n+4 )
530  nwl = iwork( 2 )
531  wu = work( n+3 )
532  wul = work( n+1 )
533  nwu = iwork( 3 )
534  END IF
535 * On exit, the interval [WL, WLU] contains a value with negcount NWL,
536 * and [WUL, WU] contains a value with negcount NWU.
537  IF( nwl.LT.0 .OR. nwl.GE.n .OR. nwu.LT.1 .OR. nwu.GT.n ) THEN
538  info = 4
539  RETURN
540  END IF
541 
542  ELSEIF( irange.EQ.valrng ) THEN
543  wl = vl
544  wu = vu
545 
546  ELSEIF( irange.EQ.allrng ) THEN
547  wl = gl
548  wu = gu
549  ENDIF
550 
551 
552 
553 * Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU.
554 * NWL accumulates the number of eigenvalues .le. WL,
555 * NWU accumulates the number of eigenvalues .le. WU
556  m = 0
557  iend = 0
558  info = 0
559  nwl = 0
560  nwu = 0
561 *
562  DO 70 jblk = 1, nsplit
563  ioff = iend
564  ibegin = ioff + 1
565  iend = isplit( jblk )
566  in = iend - ioff
567 *
568  IF( in.EQ.1 ) THEN
569 * 1x1 block
570  IF( wl.GE.d( ibegin )-pivmin )
571  $ nwl = nwl + 1
572  IF( wu.GE.d( ibegin )-pivmin )
573  $ nwu = nwu + 1
574  IF( irange.EQ.allrng .OR.
575  $ ( wl.LT.d( ibegin )-pivmin
576  $ .AND. wu.GE. d( ibegin )-pivmin ) ) THEN
577  m = m + 1
578  w( m ) = d( ibegin )
579  werr(m) = zero
580 * The gap for a single block doesn't matter for the later
581 * algorithm and is assigned an arbitrary large value
582  iblock( m ) = jblk
583  indexw( m ) = 1
584  END IF
585 
586 * Disabled 2x2 case because of a failure on the following matrix
587 * RANGE = 'I', IL = IU = 4
588 * Original Tridiagonal, d = [
589 * -0.150102010615740E+00
590 * -0.849897989384260E+00
591 * -0.128208148052635E-15
592 * 0.128257718286320E-15
593 * ];
594 * e = [
595 * -0.357171383266986E+00
596 * -0.180411241501588E-15
597 * -0.175152352710251E-15
598 * ];
599 *
600 * ELSE IF( IN.EQ.2 ) THEN
601 ** 2x2 block
602 * DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 )
603 * TMP1 = HALF*(D(IBEGIN)+D(IEND))
604 * L1 = TMP1 - DISC
605 * IF( WL.GE. L1-PIVMIN )
606 * $ NWL = NWL + 1
607 * IF( WU.GE. L1-PIVMIN )
608 * $ NWU = NWU + 1
609 * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE.
610 * $ L1-PIVMIN ) ) THEN
611 * M = M + 1
612 * W( M ) = L1
613 ** The uncertainty of eigenvalues of a 2x2 matrix is very small
614 * WERR( M ) = EPS * ABS( W( M ) ) * TWO
615 * IBLOCK( M ) = JBLK
616 * INDEXW( M ) = 1
617 * ENDIF
618 * L2 = TMP1 + DISC
619 * IF( WL.GE. L2-PIVMIN )
620 * $ NWL = NWL + 1
621 * IF( WU.GE. L2-PIVMIN )
622 * $ NWU = NWU + 1
623 * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE.
624 * $ L2-PIVMIN ) ) THEN
625 * M = M + 1
626 * W( M ) = L2
627 ** The uncertainty of eigenvalues of a 2x2 matrix is very small
628 * WERR( M ) = EPS * ABS( W( M ) ) * TWO
629 * IBLOCK( M ) = JBLK
630 * INDEXW( M ) = 2
631 * ENDIF
632  ELSE
633 * General Case - block of size IN >= 2
634 * Compute local Gerschgorin interval and use it as the initial
635 * interval for SLAEBZ
636  gu = d( ibegin )
637  gl = d( ibegin )
638  tmp1 = zero
639 
640  DO 40 j = ibegin, iend
641  gl = min( gl, gers( 2*j - 1))
642  gu = max( gu, gers(2*j) )
643  40 CONTINUE
644 * [JAN/28/2009]
645 * change SPDIAM by TNORM in lines 2 and 3 thereafter
646 * line 1: remove computation of SPDIAM (not useful anymore)
647 * SPDIAM = GU - GL
648 * GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN
649 * GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN
650  gl = gl - fudge*tnorm*eps*in - fudge*pivmin
651  gu = gu + fudge*tnorm*eps*in + fudge*pivmin
652 *
653  IF( irange.GT.1 ) THEN
654  IF( gu.LT.wl ) THEN
655 * the local block contains none of the wanted eigenvalues
656  nwl = nwl + in
657  nwu = nwu + in
658  GO TO 70
659  END IF
660 * refine search interval if possible, only range (WL,WU] matters
661  gl = max( gl, wl )
662  gu = min( gu, wu )
663  IF( gl.GE.gu )
664  $ GO TO 70
665  END IF
666 
667 * Find negcount of initial interval boundaries GL and GU
668  work( n+1 ) = gl
669  work( n+in+1 ) = gu
670  CALL slaebz( 1, 0, in, in, 1, nb, atoli, rtoli, pivmin,
671  $ d( ibegin ), e( ibegin ), e2( ibegin ),
672  $ idumma, work( n+1 ), work( n+2*in+1 ), im,
673  $ iwork, w( m+1 ), iblock( m+1 ), iinfo )
674  IF( iinfo .NE. 0 ) THEN
675  info = iinfo
676  RETURN
677  END IF
678 *
679  nwl = nwl + iwork( 1 )
680  nwu = nwu + iwork( in+1 )
681  iwoff = m - iwork( 1 )
682 
683 * Compute Eigenvalues
684  itmax = int( ( log( gu-gl+pivmin )-log( pivmin ) ) /
685  $ log( two ) ) + 2
686  CALL slaebz( 2, itmax, in, in, 1, nb, atoli, rtoli, pivmin,
687  $ d( ibegin ), e( ibegin ), e2( ibegin ),
688  $ idumma, work( n+1 ), work( n+2*in+1 ), iout,
689  $ iwork, w( m+1 ), iblock( m+1 ), iinfo )
690  IF( iinfo .NE. 0 ) THEN
691  info = iinfo
692  RETURN
693  END IF
694 *
695 * Copy eigenvalues into W and IBLOCK
696 * Use -JBLK for block number for unconverged eigenvalues.
697 * Loop over the number of output intervals from SLAEBZ
698  DO 60 j = 1, iout
699 * eigenvalue approximation is middle point of interval
700  tmp1 = half*( work( j+n )+work( j+in+n ) )
701 * semi length of error interval
702  tmp2 = half*abs( work( j+n )-work( j+in+n ) )
703  IF( j.GT.iout-iinfo ) THEN
704 * Flag non-convergence.
705  ncnvrg = .true.
706  ib = -jblk
707  ELSE
708  ib = jblk
709  END IF
710  DO 50 je = iwork( j ) + 1 + iwoff,
711  $ iwork( j+in ) + iwoff
712  w( je ) = tmp1
713  werr( je ) = tmp2
714  indexw( je ) = je - iwoff
715  iblock( je ) = ib
716  50 CONTINUE
717  60 CONTINUE
718 *
719  m = m + im
720  END IF
721  70 CONTINUE
722 
723 * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
724 * If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
725  IF( irange.EQ.indrng ) THEN
726  idiscl = il - 1 - nwl
727  idiscu = nwu - iu
728 *
729  IF( idiscl.GT.0 ) THEN
730  im = 0
731  DO 80 je = 1, m
732 * Remove some of the smallest eigenvalues from the left so that
733 * at the end IDISCL =0. Move all eigenvalues up to the left.
734  IF( w( je ).LE.wlu .AND. idiscl.GT.0 ) THEN
735  idiscl = idiscl - 1
736  ELSE
737  im = im + 1
738  w( im ) = w( je )
739  werr( im ) = werr( je )
740  indexw( im ) = indexw( je )
741  iblock( im ) = iblock( je )
742  END IF
743  80 CONTINUE
744  m = im
745  END IF
746  IF( idiscu.GT.0 ) THEN
747 * Remove some of the largest eigenvalues from the right so that
748 * at the end IDISCU =0. Move all eigenvalues up to the left.
749  im=m+1
750  DO 81 je = m, 1, -1
751  IF( w( je ).GE.wul .AND. idiscu.GT.0 ) THEN
752  idiscu = idiscu - 1
753  ELSE
754  im = im - 1
755  w( im ) = w( je )
756  werr( im ) = werr( je )
757  indexw( im ) = indexw( je )
758  iblock( im ) = iblock( je )
759  END IF
760  81 CONTINUE
761  jee = 0
762  DO 82 je = im, m
763  jee = jee + 1
764  w( jee ) = w( je )
765  werr( jee ) = werr( je )
766  indexw( jee ) = indexw( je )
767  iblock( jee ) = iblock( je )
768  82 CONTINUE
769  m = m-im+1
770  END IF
771 
772  IF( idiscl.GT.0 .OR. idiscu.GT.0 ) THEN
773 * Code to deal with effects of bad arithmetic. (If N(w) is
774 * monotone non-decreasing, this should never happen.)
775 * Some low eigenvalues to be discarded are not in (WL,WLU],
776 * or high eigenvalues to be discarded are not in (WUL,WU]
777 * so just kill off the smallest IDISCL/largest IDISCU
778 * eigenvalues, by marking the corresponding IBLOCK = 0
779  IF( idiscl.GT.0 ) THEN
780  wkill = wu
781  DO 100 jdisc = 1, idiscl
782  iw = 0
783  DO 90 je = 1, m
784  IF( iblock( je ).NE.0 .AND.
785  $ ( w( je ).LT.wkill .OR. iw.EQ.0 ) ) THEN
786  iw = je
787  wkill = w( je )
788  END IF
789  90 CONTINUE
790  iblock( iw ) = 0
791  100 CONTINUE
792  END IF
793  IF( idiscu.GT.0 ) THEN
794  wkill = wl
795  DO 120 jdisc = 1, idiscu
796  iw = 0
797  DO 110 je = 1, m
798  IF( iblock( je ).NE.0 .AND.
799  $ ( w( je ).GE.wkill .OR. iw.EQ.0 ) ) THEN
800  iw = je
801  wkill = w( je )
802  END IF
803  110 CONTINUE
804  iblock( iw ) = 0
805  120 CONTINUE
806  END IF
807 * Now erase all eigenvalues with IBLOCK set to zero
808  im = 0
809  DO 130 je = 1, m
810  IF( iblock( je ).NE.0 ) THEN
811  im = im + 1
812  w( im ) = w( je )
813  werr( im ) = werr( je )
814  indexw( im ) = indexw( je )
815  iblock( im ) = iblock( je )
816  END IF
817  130 CONTINUE
818  m = im
819  END IF
820  IF( idiscl.LT.0 .OR. idiscu.LT.0 ) THEN
821  toofew = .true.
822  END IF
823  END IF
824 *
825  IF(( irange.EQ.allrng .AND. m.NE.n ).OR.
826  $ ( irange.EQ.indrng .AND. m.NE.iu-il+1 ) ) THEN
827  toofew = .true.
828  END IF
829 
830 * If ORDER='B', do nothing the eigenvalues are already sorted by
831 * block.
832 * If ORDER='E', sort the eigenvalues from smallest to largest
833 
834  IF( lsame(order,'E') .AND. nsplit.GT.1 ) THEN
835  DO 150 je = 1, m - 1
836  ie = 0
837  tmp1 = w( je )
838  DO 140 j = je + 1, m
839  IF( w( j ).LT.tmp1 ) THEN
840  ie = j
841  tmp1 = w( j )
842  END IF
843  140 CONTINUE
844  IF( ie.NE.0 ) THEN
845  tmp2 = werr( ie )
846  itmp1 = iblock( ie )
847  itmp2 = indexw( ie )
848  w( ie ) = w( je )
849  werr( ie ) = werr( je )
850  iblock( ie ) = iblock( je )
851  indexw( ie ) = indexw( je )
852  w( je ) = tmp1
853  werr( je ) = tmp2
854  iblock( je ) = itmp1
855  indexw( je ) = itmp2
856  END IF
857  150 CONTINUE
858  END IF
859 *
860  info = 0
861  IF( ncnvrg )
862  $ info = info + 1
863  IF( toofew )
864  $ info = info + 2
865  RETURN
866 *
867 * End of SLARRD
868 *
869  END
slarrd
subroutine slarrd(RANGE, ORDER, N, VL, VU, IL, IU, GERS, RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, M, W, WERR, WL, WU, IBLOCK, INDEXW, WORK, IWORK, INFO)
SLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy.
Definition: slarrd.f:331
slaebz
subroutine slaebz(IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, NAB, WORK, IWORK, INFO)
SLAEBZ computes the number of eigenvalues of a real symmetric tridiagonal matrix which are less than ...
Definition: slaebz.f:321