LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
dlarrd.f
Go to the documentation of this file.
1 *> \brief \b DLARRD 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 DLARRD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLARRD( 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 * DOUBLE PRECISION PIVMIN, RELTOL, VL, VU, WL, WU
30 * ..
31 * .. Array Arguments ..
32 * INTEGER IBLOCK( * ), INDEXW( * ),
33 * $ ISPLIT( * ), IWORK( * )
34 * DOUBLE PRECISION D( * ), E( * ), E2( * ),
35 * $ GERS( * ), W( * ), WERR( * ), WORK( * )
36 * ..
37 *
38 *
39 *> \par Purpose:
40 * =============
41 *>
42 *> \verbatim
43 *>
44 *> DLARRD computes the eigenvalues of a symmetric tridiagonal
45 *> matrix T to suitable accuracy. This is an auxiliary code to be
46 *> called from DSTEMR.
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 DOUBLE PRECISION
95 *> \endverbatim
96 *>
97 *> \param[in] VU
98 *> \verbatim
99 *> VU is DOUBLE PRECISION
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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 DOUBLE PRECISION array, dimension (N)
191 *> On exit, the first M elements of W will contain the
192 *> eigenvalue approximations. DLARRD 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 DOUBLE PRECISION 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 DOUBLE PRECISION
209 *> \endverbatim
210 *>
211 *> \param[out] WU
212 *> \verbatim
213 *> WU is DOUBLE PRECISION
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 DLAEBZ 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. (DLARRD 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 DOUBLE PRECISION 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 DOUBLE PRECISION, 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 dlarrd( 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  DOUBLE PRECISION PIVMIN, RELTOL, VL, VU, WL, WU
333 * ..
334 * .. Array Arguments ..
335  INTEGER IBLOCK( * ), INDEXW( * ),
336  $ isplit( * ), iwork( * )
337  DOUBLE PRECISION D( * ), E( * ), E2( * ),
338  $ gers( * ), w( * ), werr( * ), work( * )
339 * ..
340 *
341 * =====================================================================
342 *
343 * .. Parameters ..
344  DOUBLE PRECISION ZERO, ONE, TWO, HALF, FUDGE
345  parameter( zero = 0.0d0, one = 1.0d0,
346  $ two = 2.0d0, 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  DOUBLE PRECISION 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  DOUBLE PRECISION DLAMCH
368  EXTERNAL lsame, ilaenv, dlamch
369 * ..
370 * .. External Subroutines ..
371  EXTERNAL dlaebz
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 = dlamch( 'P' )
429  uflow = dlamch( '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, 'DSTEBZ', ' ', 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 DLAEBZ:
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 DLAEBZ.
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 dlaebz( 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 DLAEBZ
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 dlaebz( 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 dlaebz( 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 DLAEBZ
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 DLARRD
854 *
855  END
subroutine dlaebz(IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL, RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT, NAB, WORK, IWORK, INFO)
DLAEBZ computes the number of eigenvalues of a real symmetric tridiagonal matrix which are less than ...
Definition: dlaebz.f:321
subroutine dlarrd(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)
DLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy.
Definition: dlarrd.f:323