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