LAPACK  3.6.1
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 *> 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 DOUBLE PRECISION
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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 DOUBLE PRECISION array, dimension (N)
199 *> On exit, the first M elements of W will contain the
200 *> eigenvalue approximations. DLARRD 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 DOUBLE PRECISION 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 DOUBLE PRECISION
217 *> \endverbatim
218 *>
219 *> \param[out] WU
220 *> \verbatim
221 *> WU is DOUBLE PRECISION
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 DLAEBZ 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. (DLARRD 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 DOUBLE PRECISION 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 DOUBLE PRECISION, 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 auxOTHERauxiliary
325 *
326 * =====================================================================
327  SUBROUTINE dlarrd( 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.6.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  DOUBLE PRECISION PIVMIN, RELTOL, VL, VU, WL, WU
341 * ..
342 * .. Array Arguments ..
343  INTEGER IBLOCK( * ), INDEXW( * ),
344  $ isplit( * ), iwork( * )
345  DOUBLE PRECISION D( * ), E( * ), E2( * ),
346  $ gers( * ), w( * ), werr( * ), work( * )
347 * ..
348 *
349 * =====================================================================
350 *
351 * .. Parameters ..
352  DOUBLE PRECISION ZERO, ONE, TWO, HALF, FUDGE
353  parameter ( zero = 0.0d0, one = 1.0d0,
354  $ two = 2.0d0, 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  DOUBLE PRECISION 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  DOUBLE PRECISION DLAMCH
376  EXTERNAL lsame, ilaenv, dlamch
377 * ..
378 * .. External Subroutines ..
379  EXTERNAL dlaebz
380 * ..
381 * .. Intrinsic Functions ..
382  INTRINSIC abs, int, log, max, min
383 * ..
384 * .. Executable Statements ..
385 *
386  info = 0
387 *
388 * Decode RANGE
389 *
390  IF( lsame( range, 'A' ) ) THEN
391  irange = allrng
392  ELSE IF( lsame( range, 'V' ) ) THEN
393  irange = valrng
394  ELSE IF( lsame( range, 'I' ) ) THEN
395  irange = indrng
396  ELSE
397  irange = 0
398  END IF
399 *
400 * Check for Errors
401 *
402  IF( irange.LE.0 ) THEN
403  info = -1
404  ELSE IF( .NOT.(lsame(order,'B').OR.lsame(order,'E')) ) THEN
405  info = -2
406  ELSE IF( n.LT.0 ) THEN
407  info = -3
408  ELSE IF( irange.EQ.valrng ) THEN
409  IF( vl.GE.vu )
410  $ info = -5
411  ELSE IF( irange.EQ.indrng .AND.
412  $ ( il.LT.1 .OR. il.GT.max( 1, n ) ) ) THEN
413  info = -6
414  ELSE IF( irange.EQ.indrng .AND.
415  $ ( iu.LT.min( n, il ) .OR. iu.GT.n ) ) THEN
416  info = -7
417  END IF
418 *
419  IF( info.NE.0 ) THEN
420  RETURN
421  END IF
422 
423 * Initialize error flags
424  info = 0
425  ncnvrg = .false.
426  toofew = .false.
427 
428 * Quick return if possible
429  m = 0
430  IF( n.EQ.0 ) RETURN
431 
432 * Simplification:
433  IF( irange.EQ.indrng .AND. il.EQ.1 .AND. iu.EQ.n ) irange = 1
434 
435 * Get machine constants
436  eps = dlamch( 'P' )
437  uflow = dlamch( 'U' )
438 
439 
440 * Special Case when N=1
441 * Treat case of 1x1 matrix for quick return
442  IF( n.EQ.1 ) THEN
443  IF( (irange.EQ.allrng).OR.
444  $ ((irange.EQ.valrng).AND.(d(1).GT.vl).AND.(d(1).LE.vu)).OR.
445  $ ((irange.EQ.indrng).AND.(il.EQ.1).AND.(iu.EQ.1)) ) THEN
446  m = 1
447  w(1) = d(1)
448 * The computation error of the eigenvalue is zero
449  werr(1) = zero
450  iblock( 1 ) = 1
451  indexw( 1 ) = 1
452  ENDIF
453  RETURN
454  END IF
455 
456 * NB is the minimum vector length for vector bisection, or 0
457 * if only scalar is to be done.
458  nb = ilaenv( 1, 'DSTEBZ', ' ', n, -1, -1, -1 )
459  IF( nb.LE.1 ) nb = 0
460 
461 * Find global spectral radius
462  gl = d(1)
463  gu = d(1)
464  DO 5 i = 1,n
465  gl = min( gl, gers( 2*i - 1))
466  gu = max( gu, gers(2*i) )
467  5 CONTINUE
468 * Compute global Gerschgorin bounds and spectral diameter
469  tnorm = max( abs( gl ), abs( gu ) )
470  gl = gl - fudge*tnorm*eps*n - fudge*two*pivmin
471  gu = gu + fudge*tnorm*eps*n + fudge*two*pivmin
472 * [JAN/28/2009] remove the line below since SPDIAM variable not use
473 * SPDIAM = GU - GL
474 * Input arguments for DLAEBZ:
475 * The relative tolerance. An interval (a,b] lies within
476 * "relative tolerance" if b-a < RELTOL*max(|a|,|b|),
477  rtoli = reltol
478 * Set the absolute tolerance for interval convergence to zero to force
479 * interval convergence based on relative size of the interval.
480 * This is dangerous because intervals might not converge when RELTOL is
481 * small. But at least a very small number should be selected so that for
482 * strongly graded matrices, the code can get relatively accurate
483 * eigenvalues.
484  atoli = fudge*two*uflow + fudge*two*pivmin
485 
486  IF( irange.EQ.indrng ) THEN
487 
488 * RANGE='I': Compute an interval containing eigenvalues
489 * IL through IU. The initial interval [GL,GU] from the global
490 * Gerschgorin bounds GL and GU is refined by DLAEBZ.
491  itmax = int( ( log( tnorm+pivmin )-log( pivmin ) ) /
492  $ log( two ) ) + 2
493  work( n+1 ) = gl
494  work( n+2 ) = gl
495  work( n+3 ) = gu
496  work( n+4 ) = gu
497  work( n+5 ) = gl
498  work( n+6 ) = gu
499  iwork( 1 ) = -1
500  iwork( 2 ) = -1
501  iwork( 3 ) = n + 1
502  iwork( 4 ) = n + 1
503  iwork( 5 ) = il - 1
504  iwork( 6 ) = iu
505 *
506  CALL dlaebz( 3, itmax, n, 2, 2, nb, atoli, rtoli, pivmin,
507  $ d, e, e2, iwork( 5 ), work( n+1 ), work( n+5 ), iout,
508  $ iwork, w, iblock, iinfo )
509  IF( iinfo .NE. 0 ) THEN
510  info = iinfo
511  RETURN
512  END IF
513 * On exit, output intervals may not be ordered by ascending negcount
514  IF( iwork( 6 ).EQ.iu ) THEN
515  wl = work( n+1 )
516  wlu = work( n+3 )
517  nwl = iwork( 1 )
518  wu = work( n+4 )
519  wul = work( n+2 )
520  nwu = iwork( 4 )
521  ELSE
522  wl = work( n+2 )
523  wlu = work( n+4 )
524  nwl = iwork( 2 )
525  wu = work( n+3 )
526  wul = work( n+1 )
527  nwu = iwork( 3 )
528  END IF
529 * On exit, the interval [WL, WLU] contains a value with negcount NWL,
530 * and [WUL, WU] contains a value with negcount NWU.
531  IF( nwl.LT.0 .OR. nwl.GE.n .OR. nwu.LT.1 .OR. nwu.GT.n ) THEN
532  info = 4
533  RETURN
534  END IF
535 
536  ELSEIF( irange.EQ.valrng ) THEN
537  wl = vl
538  wu = vu
539 
540  ELSEIF( irange.EQ.allrng ) THEN
541  wl = gl
542  wu = gu
543  ENDIF
544 
545 
546 
547 * Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU.
548 * NWL accumulates the number of eigenvalues .le. WL,
549 * NWU accumulates the number of eigenvalues .le. WU
550  m = 0
551  iend = 0
552  info = 0
553  nwl = 0
554  nwu = 0
555 *
556  DO 70 jblk = 1, nsplit
557  ioff = iend
558  ibegin = ioff + 1
559  iend = isplit( jblk )
560  in = iend - ioff
561 *
562  IF( in.EQ.1 ) THEN
563 * 1x1 block
564  IF( wl.GE.d( ibegin )-pivmin )
565  $ nwl = nwl + 1
566  IF( wu.GE.d( ibegin )-pivmin )
567  $ nwu = nwu + 1
568  IF( irange.EQ.allrng .OR.
569  $ ( wl.LT.d( ibegin )-pivmin
570  $ .AND. wu.GE. d( ibegin )-pivmin ) ) THEN
571  m = m + 1
572  w( m ) = d( ibegin )
573  werr(m) = zero
574 * The gap for a single block doesn't matter for the later
575 * algorithm and is assigned an arbitrary large value
576  iblock( m ) = jblk
577  indexw( m ) = 1
578  END IF
579 
580 * Disabled 2x2 case because of a failure on the following matrix
581 * RANGE = 'I', IL = IU = 4
582 * Original Tridiagonal, d = [
583 * -0.150102010615740E+00
584 * -0.849897989384260E+00
585 * -0.128208148052635E-15
586 * 0.128257718286320E-15
587 * ];
588 * e = [
589 * -0.357171383266986E+00
590 * -0.180411241501588E-15
591 * -0.175152352710251E-15
592 * ];
593 *
594 * ELSE IF( IN.EQ.2 ) THEN
595 ** 2x2 block
596 * DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 )
597 * TMP1 = HALF*(D(IBEGIN)+D(IEND))
598 * L1 = TMP1 - DISC
599 * IF( WL.GE. L1-PIVMIN )
600 * $ NWL = NWL + 1
601 * IF( WU.GE. L1-PIVMIN )
602 * $ NWU = NWU + 1
603 * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE.
604 * $ L1-PIVMIN ) ) THEN
605 * M = M + 1
606 * W( M ) = L1
607 ** The uncertainty of eigenvalues of a 2x2 matrix is very small
608 * WERR( M ) = EPS * ABS( W( M ) ) * TWO
609 * IBLOCK( M ) = JBLK
610 * INDEXW( M ) = 1
611 * ENDIF
612 * L2 = TMP1 + DISC
613 * IF( WL.GE. L2-PIVMIN )
614 * $ NWL = NWL + 1
615 * IF( WU.GE. L2-PIVMIN )
616 * $ NWU = NWU + 1
617 * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE.
618 * $ L2-PIVMIN ) ) THEN
619 * M = M + 1
620 * W( M ) = L2
621 ** The uncertainty of eigenvalues of a 2x2 matrix is very small
622 * WERR( M ) = EPS * ABS( W( M ) ) * TWO
623 * IBLOCK( M ) = JBLK
624 * INDEXW( M ) = 2
625 * ENDIF
626  ELSE
627 * General Case - block of size IN >= 2
628 * Compute local Gerschgorin interval and use it as the initial
629 * interval for DLAEBZ
630  gu = d( ibegin )
631  gl = d( ibegin )
632  tmp1 = zero
633 
634  DO 40 j = ibegin, iend
635  gl = min( gl, gers( 2*j - 1))
636  gu = max( gu, gers(2*j) )
637  40 CONTINUE
638 * [JAN/28/2009]
639 * change SPDIAM by TNORM in lines 2 and 3 thereafter
640 * line 1: remove computation of SPDIAM (not useful anymore)
641 * SPDIAM = GU - GL
642 * GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN
643 * GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN
644  gl = gl - fudge*tnorm*eps*in - fudge*pivmin
645  gu = gu + fudge*tnorm*eps*in + fudge*pivmin
646 *
647  IF( irange.GT.1 ) THEN
648  IF( gu.LT.wl ) THEN
649 * the local block contains none of the wanted eigenvalues
650  nwl = nwl + in
651  nwu = nwu + in
652  GO TO 70
653  END IF
654 * refine search interval if possible, only range (WL,WU] matters
655  gl = max( gl, wl )
656  gu = min( gu, wu )
657  IF( gl.GE.gu )
658  $ GO TO 70
659  END IF
660 
661 * Find negcount of initial interval boundaries GL and GU
662  work( n+1 ) = gl
663  work( n+in+1 ) = gu
664  CALL dlaebz( 1, 0, in, in, 1, nb, atoli, rtoli, pivmin,
665  $ d( ibegin ), e( ibegin ), e2( ibegin ),
666  $ idumma, work( n+1 ), work( n+2*in+1 ), im,
667  $ iwork, w( m+1 ), iblock( m+1 ), iinfo )
668  IF( iinfo .NE. 0 ) THEN
669  info = iinfo
670  RETURN
671  END IF
672 *
673  nwl = nwl + iwork( 1 )
674  nwu = nwu + iwork( in+1 )
675  iwoff = m - iwork( 1 )
676 
677 * Compute Eigenvalues
678  itmax = int( ( log( gu-gl+pivmin )-log( pivmin ) ) /
679  $ log( two ) ) + 2
680  CALL dlaebz( 2, itmax, in, in, 1, nb, atoli, rtoli, pivmin,
681  $ d( ibegin ), e( ibegin ), e2( ibegin ),
682  $ idumma, work( n+1 ), work( n+2*in+1 ), iout,
683  $ iwork, w( m+1 ), iblock( m+1 ), iinfo )
684  IF( iinfo .NE. 0 ) THEN
685  info = iinfo
686  RETURN
687  END IF
688 *
689 * Copy eigenvalues into W and IBLOCK
690 * Use -JBLK for block number for unconverged eigenvalues.
691 * Loop over the number of output intervals from DLAEBZ
692  DO 60 j = 1, iout
693 * eigenvalue approximation is middle point of interval
694  tmp1 = half*( work( j+n )+work( j+in+n ) )
695 * semi length of error interval
696  tmp2 = half*abs( work( j+n )-work( j+in+n ) )
697  IF( j.GT.iout-iinfo ) THEN
698 * Flag non-convergence.
699  ncnvrg = .true.
700  ib = -jblk
701  ELSE
702  ib = jblk
703  END IF
704  DO 50 je = iwork( j ) + 1 + iwoff,
705  $ iwork( j+in ) + iwoff
706  w( je ) = tmp1
707  werr( je ) = tmp2
708  indexw( je ) = je - iwoff
709  iblock( je ) = ib
710  50 CONTINUE
711  60 CONTINUE
712 *
713  m = m + im
714  END IF
715  70 CONTINUE
716 
717 * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
718 * If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
719  IF( irange.EQ.indrng ) THEN
720  idiscl = il - 1 - nwl
721  idiscu = nwu - iu
722 *
723  IF( idiscl.GT.0 ) THEN
724  im = 0
725  DO 80 je = 1, m
726 * Remove some of the smallest eigenvalues from the left so that
727 * at the end IDISCL =0. Move all eigenvalues up to the left.
728  IF( w( je ).LE.wlu .AND. idiscl.GT.0 ) THEN
729  idiscl = idiscl - 1
730  ELSE
731  im = im + 1
732  w( im ) = w( je )
733  werr( im ) = werr( je )
734  indexw( im ) = indexw( je )
735  iblock( im ) = iblock( je )
736  END IF
737  80 CONTINUE
738  m = im
739  END IF
740  IF( idiscu.GT.0 ) THEN
741 * Remove some of the largest eigenvalues from the right so that
742 * at the end IDISCU =0. Move all eigenvalues up to the left.
743  im=m+1
744  DO 81 je = m, 1, -1
745  IF( w( je ).GE.wul .AND. idiscu.GT.0 ) THEN
746  idiscu = idiscu - 1
747  ELSE
748  im = im - 1
749  w( im ) = w( je )
750  werr( im ) = werr( je )
751  indexw( im ) = indexw( je )
752  iblock( im ) = iblock( je )
753  END IF
754  81 CONTINUE
755  jee = 0
756  DO 82 je = im, m
757  jee = jee + 1
758  w( jee ) = w( je )
759  werr( jee ) = werr( je )
760  indexw( jee ) = indexw( je )
761  iblock( jee ) = iblock( je )
762  82 CONTINUE
763  m = m-im+1
764  END IF
765 
766  IF( idiscl.GT.0 .OR. idiscu.GT.0 ) THEN
767 * Code to deal with effects of bad arithmetic. (If N(w) is
768 * monotone non-decreasing, this should never happen.)
769 * Some low eigenvalues to be discarded are not in (WL,WLU],
770 * or high eigenvalues to be discarded are not in (WUL,WU]
771 * so just kill off the smallest IDISCL/largest IDISCU
772 * eigenvalues, by marking the corresponding IBLOCK = 0
773  IF( idiscl.GT.0 ) THEN
774  wkill = wu
775  DO 100 jdisc = 1, idiscl
776  iw = 0
777  DO 90 je = 1, m
778  IF( iblock( je ).NE.0 .AND.
779  $ ( w( je ).LT.wkill .OR. iw.EQ.0 ) ) THEN
780  iw = je
781  wkill = w( je )
782  END IF
783  90 CONTINUE
784  iblock( iw ) = 0
785  100 CONTINUE
786  END IF
787  IF( idiscu.GT.0 ) THEN
788  wkill = wl
789  DO 120 jdisc = 1, idiscu
790  iw = 0
791  DO 110 je = 1, m
792  IF( iblock( je ).NE.0 .AND.
793  $ ( w( je ).GE.wkill .OR. iw.EQ.0 ) ) THEN
794  iw = je
795  wkill = w( je )
796  END IF
797  110 CONTINUE
798  iblock( iw ) = 0
799  120 CONTINUE
800  END IF
801 * Now erase all eigenvalues with IBLOCK set to zero
802  im = 0
803  DO 130 je = 1, m
804  IF( iblock( je ).NE.0 ) THEN
805  im = im + 1
806  w( im ) = w( je )
807  werr( im ) = werr( je )
808  indexw( im ) = indexw( je )
809  iblock( im ) = iblock( je )
810  END IF
811  130 CONTINUE
812  m = im
813  END IF
814  IF( idiscl.LT.0 .OR. idiscu.LT.0 ) THEN
815  toofew = .true.
816  END IF
817  END IF
818 *
819  IF(( irange.EQ.allrng .AND. m.NE.n ).OR.
820  $ ( irange.EQ.indrng .AND. m.NE.iu-il+1 ) ) THEN
821  toofew = .true.
822  END IF
823 
824 * If ORDER='B', do nothing the eigenvalues are already sorted by
825 * block.
826 * If ORDER='E', sort the eigenvalues from smallest to largest
827 
828  IF( lsame(order,'E') .AND. nsplit.GT.1 ) THEN
829  DO 150 je = 1, m - 1
830  ie = 0
831  tmp1 = w( je )
832  DO 140 j = je + 1, m
833  IF( w( j ).LT.tmp1 ) THEN
834  ie = j
835  tmp1 = w( j )
836  END IF
837  140 CONTINUE
838  IF( ie.NE.0 ) THEN
839  tmp2 = werr( ie )
840  itmp1 = iblock( ie )
841  itmp2 = indexw( ie )
842  w( ie ) = w( je )
843  werr( ie ) = werr( je )
844  iblock( ie ) = iblock( je )
845  indexw( ie ) = indexw( je )
846  w( je ) = tmp1
847  werr( je ) = tmp2
848  iblock( je ) = itmp1
849  indexw( je ) = itmp2
850  END IF
851  150 CONTINUE
852  END IF
853 *
854  info = 0
855  IF( ncnvrg )
856  $ info = info + 1
857  IF( toofew )
858  $ info = info + 2
859  RETURN
860 *
861 * End of DLARRD
862 *
863  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:331