LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
dlarre.f
Go to the documentation of this file.
1 *> \brief \b DLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduced block Ti, finds base representations and eigenvalues.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DLARRE + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarre.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarre.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarre.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2,
22 * RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
23 * W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
24 * WORK, IWORK, INFO )
25 *
26 * .. Scalar Arguments ..
27 * CHARACTER RANGE
28 * INTEGER IL, INFO, IU, M, N, NSPLIT
29 * DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
30 * ..
31 * .. Array Arguments ..
32 * INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ),
33 * $ INDEXW( * )
34 * DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ),
35 * $ W( * ),WERR( * ), WGAP( * ), WORK( * )
36 * ..
37 *
38 *
39 *> \par Purpose:
40 * =============
41 *>
42 *> \verbatim
43 *>
44 *> To find the desired eigenvalues of a given real symmetric
45 *> tridiagonal matrix T, DLARRE sets any "small" off-diagonal
46 *> elements to zero, and for each unreduced block T_i, it finds
47 *> (a) a suitable shift at one end of the block's spectrum,
48 *> (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
49 *> (c) eigenvalues of each L_i D_i L_i^T.
50 *> The representations and eigenvalues found are then used by
51 *> DSTEMR to compute the eigenvectors of T.
52 *> The accuracy varies depending on whether bisection is used to
53 *> find a few eigenvalues or the dqds algorithm (subroutine DLASQ2) to
54 *> conpute all and then discard any unwanted one.
55 *> As an added benefit, DLARRE also outputs the n
56 *> Gerschgorin intervals for the matrices L_i D_i L_i^T.
57 *> \endverbatim
58 *
59 * Arguments:
60 * ==========
61 *
62 *> \param[in] RANGE
63 *> \verbatim
64 *> RANGE is CHARACTER*1
65 *> = 'A': ("All") all eigenvalues will be found.
66 *> = 'V': ("Value") all eigenvalues in the half-open interval
67 *> (VL, VU] will be found.
68 *> = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
69 *> entire matrix) will be found.
70 *> \endverbatim
71 *>
72 *> \param[in] N
73 *> \verbatim
74 *> N is INTEGER
75 *> The order of the matrix. N > 0.
76 *> \endverbatim
77 *>
78 *> \param[in,out] VL
79 *> \verbatim
80 *> VL is DOUBLE PRECISION
81 *> \endverbatim
82 *>
83 *> \param[in,out] VU
84 *> \verbatim
85 *> VU is DOUBLE PRECISION
86 *> If RANGE='V', the lower and upper bounds for the eigenvalues.
87 *> Eigenvalues less than or equal to VL, or greater than VU,
88 *> will not be returned. VL < VU.
89 *> If RANGE='I' or ='A', DLARRE computes bounds on the desired
90 *> part of the spectrum.
91 *> \endverbatim
92 *>
93 *> \param[in] IL
94 *> \verbatim
95 *> IL is INTEGER
96 *> \endverbatim
97 *>
98 *> \param[in] IU
99 *> \verbatim
100 *> IU is INTEGER
101 *> If RANGE='I', the indices (in ascending order) of the
102 *> smallest and largest eigenvalues to be returned.
103 *> 1 <= IL <= IU <= N.
104 *> \endverbatim
105 *>
106 *> \param[in,out] D
107 *> \verbatim
108 *> D is DOUBLE PRECISION array, dimension (N)
109 *> On entry, the N diagonal elements of the tridiagonal
110 *> matrix T.
111 *> On exit, the N diagonal elements of the diagonal
112 *> matrices D_i.
113 *> \endverbatim
114 *>
115 *> \param[in,out] E
116 *> \verbatim
117 *> E is DOUBLE PRECISION array, dimension (N)
118 *> On entry, the first (N-1) entries contain the subdiagonal
119 *> elements of the tridiagonal matrix T; E(N) need not be set.
120 *> On exit, E contains the subdiagonal elements of the unit
121 *> bidiagonal matrices L_i. The entries E( ISPLIT( I ) ),
122 *> 1 <= I <= NSPLIT, contain the base points sigma_i on output.
123 *> \endverbatim
124 *>
125 *> \param[in,out] E2
126 *> \verbatim
127 *> E2 is DOUBLE PRECISION array, dimension (N)
128 *> On entry, the first (N-1) entries contain the SQUARES of the
129 *> subdiagonal elements of the tridiagonal matrix T;
130 *> E2(N) need not be set.
131 *> On exit, the entries E2( ISPLIT( I ) ),
132 *> 1 <= I <= NSPLIT, have been set to zero
133 *> \endverbatim
134 *>
135 *> \param[in] RTOL1
136 *> \verbatim
137 *> RTOL1 is DOUBLE PRECISION
138 *> \endverbatim
139 *>
140 *> \param[in] RTOL2
141 *> \verbatim
142 *> RTOL2 is DOUBLE PRECISION
143 *> Parameters for bisection.
144 *> An interval [LEFT,RIGHT] has converged if
145 *> RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
146 *> \endverbatim
147 *>
148 *> \param[in] SPLTOL
149 *> \verbatim
150 *> SPLTOL is DOUBLE PRECISION
151 *> The threshold for splitting.
152 *> \endverbatim
153 *>
154 *> \param[out] NSPLIT
155 *> \verbatim
156 *> NSPLIT is INTEGER
157 *> The number of blocks T splits into. 1 <= NSPLIT <= N.
158 *> \endverbatim
159 *>
160 *> \param[out] ISPLIT
161 *> \verbatim
162 *> ISPLIT is INTEGER array, dimension (N)
163 *> The splitting points, at which T breaks up into blocks.
164 *> The first block consists of rows/columns 1 to ISPLIT(1),
165 *> the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
166 *> etc., and the NSPLIT-th consists of rows/columns
167 *> ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
168 *> \endverbatim
169 *>
170 *> \param[out] M
171 *> \verbatim
172 *> M is INTEGER
173 *> The total number of eigenvalues (of all L_i D_i L_i^T)
174 *> found.
175 *> \endverbatim
176 *>
177 *> \param[out] W
178 *> \verbatim
179 *> W is DOUBLE PRECISION array, dimension (N)
180 *> The first M elements contain the eigenvalues. The
181 *> eigenvalues of each of the blocks, L_i D_i L_i^T, are
182 *> sorted in ascending order ( DLARRE may use the
183 *> remaining N-M elements as workspace).
184 *> \endverbatim
185 *>
186 *> \param[out] WERR
187 *> \verbatim
188 *> WERR is DOUBLE PRECISION array, dimension (N)
189 *> The error bound on the corresponding eigenvalue in W.
190 *> \endverbatim
191 *>
192 *> \param[out] WGAP
193 *> \verbatim
194 *> WGAP is DOUBLE PRECISION array, dimension (N)
195 *> The separation from the right neighbor eigenvalue in W.
196 *> The gap is only with respect to the eigenvalues of the same block
197 *> as each block has its own representation tree.
198 *> Exception: at the right end of a block we store the left gap
199 *> \endverbatim
200 *>
201 *> \param[out] IBLOCK
202 *> \verbatim
203 *> IBLOCK is INTEGER array, dimension (N)
204 *> The indices of the blocks (submatrices) associated with the
205 *> corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
206 *> W(i) belongs to the first block from the top, =2 if W(i)
207 *> belongs to the second block, etc.
208 *> \endverbatim
209 *>
210 *> \param[out] INDEXW
211 *> \verbatim
212 *> INDEXW is INTEGER array, dimension (N)
213 *> The indices of the eigenvalues within each block (submatrix);
214 *> for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
215 *> i-th eigenvalue W(i) is the 10-th eigenvalue in block 2
216 *> \endverbatim
217 *>
218 *> \param[out] GERS
219 *> \verbatim
220 *> GERS is DOUBLE PRECISION array, dimension (2*N)
221 *> The N Gerschgorin intervals (the i-th Gerschgorin interval
222 *> is (GERS(2*i-1), GERS(2*i)).
223 *> \endverbatim
224 *>
225 *> \param[out] PIVMIN
226 *> \verbatim
227 *> PIVMIN is DOUBLE PRECISION
228 *> The minimum pivot in the Sturm sequence for T.
229 *> \endverbatim
230 *>
231 *> \param[out] WORK
232 *> \verbatim
233 *> WORK is DOUBLE PRECISION array, dimension (6*N)
234 *> Workspace.
235 *> \endverbatim
236 *>
237 *> \param[out] IWORK
238 *> \verbatim
239 *> IWORK is INTEGER array, dimension (5*N)
240 *> Workspace.
241 *> \endverbatim
242 *>
243 *> \param[out] INFO
244 *> \verbatim
245 *> INFO is INTEGER
246 *> = 0: successful exit
247 *> > 0: A problem occured in DLARRE.
248 *> < 0: One of the called subroutines signaled an internal problem.
249 *> Needs inspection of the corresponding parameter IINFO
250 *> for further information.
251 *>
252 *> =-1: Problem in DLARRD.
253 *> = 2: No base representation could be found in MAXTRY iterations.
254 *> Increasing MAXTRY and recompilation might be a remedy.
255 *> =-3: Problem in DLARRB when computing the refined root
256 *> representation for DLASQ2.
257 *> =-4: Problem in DLARRB when preforming bisection on the
258 *> desired part of the spectrum.
259 *> =-5: Problem in DLASQ2.
260 *> =-6: Problem in DLASQ2.
261 *> \endverbatim
262 *
263 * Authors:
264 * ========
265 *
266 *> \author Univ. of Tennessee
267 *> \author Univ. of California Berkeley
268 *> \author Univ. of Colorado Denver
269 *> \author NAG Ltd.
270 *
271 *> \date September 2012
272 *
273 *> \ingroup auxOTHERauxiliary
274 *
275 *> \par Further Details:
276 * =====================
277 *>
278 *> \verbatim
279 *>
280 *> The base representations are required to suffer very little
281 *> element growth and consequently define all their eigenvalues to
282 *> high relative accuracy.
283 *> \endverbatim
284 *
285 *> \par Contributors:
286 * ==================
287 *>
288 *> Beresford Parlett, University of California, Berkeley, USA \n
289 *> Jim Demmel, University of California, Berkeley, USA \n
290 *> Inderjit Dhillon, University of Texas, Austin, USA \n
291 *> Osni Marques, LBNL/NERSC, USA \n
292 *> Christof Voemel, University of California, Berkeley, USA \n
293 *>
294 * =====================================================================
295  SUBROUTINE dlarre( RANGE, N, VL, VU, IL, IU, D, E, E2,
296  $ rtol1, rtol2, spltol, nsplit, isplit, m,
297  $ w, werr, wgap, iblock, indexw, gers, pivmin,
298  $ work, iwork, info )
299 *
300 * -- LAPACK auxiliary routine (version 3.4.2) --
301 * -- LAPACK is a software package provided by Univ. of Tennessee, --
302 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
303 * September 2012
304 *
305 * .. Scalar Arguments ..
306  CHARACTER RANGE
307  INTEGER IL, INFO, IU, M, N, NSPLIT
308  DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
309 * ..
310 * .. Array Arguments ..
311  INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ),
312  $ indexw( * )
313  DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ),
314  $ w( * ),werr( * ), wgap( * ), work( * )
315 * ..
316 *
317 * =====================================================================
318 *
319 * .. Parameters ..
320  DOUBLE PRECISION FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
321  $ maxgrowth, one, pert, two, zero
322  parameter( zero = 0.0d0, one = 1.0d0,
323  $ two = 2.0d0, four=4.0d0,
324  $ hndrd = 100.0d0,
325  $ pert = 8.0d0,
326  $ half = one/two, fourth = one/four, fac= half,
327  $ maxgrowth = 64.0d0, fudge = 2.0d0 )
328  INTEGER MAXTRY, ALLRNG, INDRNG, VALRNG
329  parameter( maxtry = 6, allrng = 1, indrng = 2,
330  $ valrng = 3 )
331 * ..
332 * .. Local Scalars ..
333  LOGICAL FORCEB, NOREP, USEDQD
334  INTEGER CNT, CNT1, CNT2, I, IBEGIN, IDUM, IEND, IINFO,
335  $ in, indl, indu, irange, j, jblk, mb, mm,
336  $ wbegin, wend
337  DOUBLE PRECISION AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS,
338  $ emax, eold, eps, gl, gu, isleft, isrght, rtl,
339  $ rtol, s1, s2, safmin, sgndef, sigma, spdiam,
340  $ tau, tmp, tmp1
341 
342 
343 * ..
344 * .. Local Arrays ..
345  INTEGER ISEED( 4 )
346 * ..
347 * .. External Functions ..
348  LOGICAL LSAME
349  DOUBLE PRECISION DLAMCH
350  EXTERNAL dlamch, lsame
351 
352 * ..
353 * .. External Subroutines ..
354  EXTERNAL dcopy, dlarnv, dlarra, dlarrb, dlarrc, dlarrd,
355  $ dlasq2
356 * ..
357 * .. Intrinsic Functions ..
358  INTRINSIC abs, max, min
359 
360 * ..
361 * .. Executable Statements ..
362 *
363 
364  info = 0
365 
366 *
367 * Decode RANGE
368 *
369  IF( lsame( range, 'A' ) ) THEN
370  irange = allrng
371  ELSE IF( lsame( range, 'V' ) ) THEN
372  irange = valrng
373  ELSE IF( lsame( range, 'I' ) ) THEN
374  irange = indrng
375  END IF
376 
377  m = 0
378 
379 * Get machine constants
380  safmin = dlamch( 'S' )
381  eps = dlamch( 'P' )
382 
383 * Set parameters
384  rtl = sqrt(eps)
385  bsrtol = sqrt(eps)
386 
387 * Treat case of 1x1 matrix for quick return
388  IF( n.EQ.1 ) THEN
389  IF( (irange.EQ.allrng).OR.
390  $ ((irange.EQ.valrng).AND.(d(1).GT.vl).AND.(d(1).LE.vu)).OR.
391  $ ((irange.EQ.indrng).AND.(il.EQ.1).AND.(iu.EQ.1)) ) THEN
392  m = 1
393  w(1) = d(1)
394 * The computation error of the eigenvalue is zero
395  werr(1) = zero
396  wgap(1) = zero
397  iblock( 1 ) = 1
398  indexw( 1 ) = 1
399  gers(1) = d( 1 )
400  gers(2) = d( 1 )
401  ENDIF
402 * store the shift for the initial RRR, which is zero in this case
403  e(1) = zero
404  RETURN
405  END IF
406 
407 * General case: tridiagonal matrix of order > 1
408 *
409 * Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter.
410 * Compute maximum off-diagonal entry and pivmin.
411  gl = d(1)
412  gu = d(1)
413  eold = zero
414  emax = zero
415  e(n) = zero
416  DO 5 i = 1,n
417  werr(i) = zero
418  wgap(i) = zero
419  eabs = abs( e(i) )
420  IF( eabs .GE. emax ) THEN
421  emax = eabs
422  END IF
423  tmp1 = eabs + eold
424  gers( 2*i-1) = d(i) - tmp1
425  gl = min( gl, gers( 2*i - 1))
426  gers( 2*i ) = d(i) + tmp1
427  gu = max( gu, gers(2*i) )
428  eold = eabs
429  5 CONTINUE
430 * The minimum pivot allowed in the Sturm sequence for T
431  pivmin = safmin * max( one, emax**2 )
432 * Compute spectral diameter. The Gerschgorin bounds give an
433 * estimate that is wrong by at most a factor of SQRT(2)
434  spdiam = gu - gl
435 
436 * Compute splitting points
437  CALL dlarra( n, d, e, e2, spltol, spdiam,
438  $ nsplit, isplit, iinfo )
439 
440 * Can force use of bisection instead of faster DQDS.
441 * Option left in the code for future multisection work.
442  forceb = .false.
443 
444 * Initialize USEDQD, DQDS should be used for ALLRNG unless someone
445 * explicitly wants bisection.
446  usedqd = (( irange.EQ.allrng ) .AND. (.NOT.forceb))
447 
448  IF( (irange.EQ.allrng) .AND. (.NOT. forceb) ) THEN
449 * Set interval [VL,VU] that contains all eigenvalues
450  vl = gl
451  vu = gu
452  ELSE
453 * We call DLARRD to find crude approximations to the eigenvalues
454 * in the desired range. In case IRANGE = INDRNG, we also obtain the
455 * interval (VL,VU] that contains all the wanted eigenvalues.
456 * An interval [LEFT,RIGHT] has converged if
457 * RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT))
458 * DLARRD needs a WORK of size 4*N, IWORK of size 3*N
459  CALL dlarrd( range, 'B', n, vl, vu, il, iu, gers,
460  $ bsrtol, d, e, e2, pivmin, nsplit, isplit,
461  $ mm, w, werr, vl, vu, iblock, indexw,
462  $ work, iwork, iinfo )
463  IF( iinfo.NE.0 ) THEN
464  info = -1
465  RETURN
466  ENDIF
467 * Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0
468  DO 14 i = mm+1,n
469  w( i ) = zero
470  werr( i ) = zero
471  iblock( i ) = 0
472  indexw( i ) = 0
473  14 CONTINUE
474  END IF
475 
476 
477 ***
478 * Loop over unreduced blocks
479  ibegin = 1
480  wbegin = 1
481  DO 170 jblk = 1, nsplit
482  iend = isplit( jblk )
483  in = iend - ibegin + 1
484 
485 * 1 X 1 block
486  IF( in.EQ.1 ) THEN
487  IF( (irange.EQ.allrng).OR.( (irange.EQ.valrng).AND.
488  $ ( d( ibegin ).GT.vl ).AND.( d( ibegin ).LE.vu ) )
489  $ .OR. ( (irange.EQ.indrng).AND.(iblock(wbegin).EQ.jblk))
490  $ ) THEN
491  m = m + 1
492  w( m ) = d( ibegin )
493  werr(m) = zero
494 * The gap for a single block doesn't matter for the later
495 * algorithm and is assigned an arbitrary large value
496  wgap(m) = zero
497  iblock( m ) = jblk
498  indexw( m ) = 1
499  wbegin = wbegin + 1
500  ENDIF
501 * E( IEND ) holds the shift for the initial RRR
502  e( iend ) = zero
503  ibegin = iend + 1
504  GO TO 170
505  END IF
506 *
507 * Blocks of size larger than 1x1
508 *
509 * E( IEND ) will hold the shift for the initial RRR, for now set it =0
510  e( iend ) = zero
511 *
512 * Find local outer bounds GL,GU for the block
513  gl = d(ibegin)
514  gu = d(ibegin)
515  DO 15 i = ibegin , iend
516  gl = min( gers( 2*i-1 ), gl )
517  gu = max( gers( 2*i ), gu )
518  15 CONTINUE
519  spdiam = gu - gl
520 
521  IF(.NOT. ((irange.EQ.allrng).AND.(.NOT.forceb)) ) THEN
522 * Count the number of eigenvalues in the current block.
523  mb = 0
524  DO 20 i = wbegin,mm
525  IF( iblock(i).EQ.jblk ) THEN
526  mb = mb+1
527  ELSE
528  GOTO 21
529  ENDIF
530  20 CONTINUE
531  21 CONTINUE
532 
533  IF( mb.EQ.0) THEN
534 * No eigenvalue in the current block lies in the desired range
535 * E( IEND ) holds the shift for the initial RRR
536  e( iend ) = zero
537  ibegin = iend + 1
538  GO TO 170
539  ELSE
540 
541 * Decide whether dqds or bisection is more efficient
542  usedqd = ( (mb .GT. fac*in) .AND. (.NOT.forceb) )
543  wend = wbegin + mb - 1
544 * Calculate gaps for the current block
545 * In later stages, when representations for individual
546 * eigenvalues are different, we use SIGMA = E( IEND ).
547  sigma = zero
548  DO 30 i = wbegin, wend - 1
549  wgap( i ) = max( zero,
550  $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
551  30 CONTINUE
552  wgap( wend ) = max( zero,
553  $ vu - sigma - (w( wend )+werr( wend )))
554 * Find local index of the first and last desired evalue.
555  indl = indexw(wbegin)
556  indu = indexw( wend )
557  ENDIF
558  ENDIF
559  IF(( (irange.EQ.allrng) .AND. (.NOT. forceb) ).OR.usedqd) THEN
560 * Case of DQDS
561 * Find approximations to the extremal eigenvalues of the block
562  CALL dlarrk( in, 1, gl, gu, d(ibegin),
563  $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
564  IF( iinfo.NE.0 ) THEN
565  info = -1
566  RETURN
567  ENDIF
568  isleft = max(gl, tmp - tmp1
569  $ - hndrd * eps* abs(tmp - tmp1))
570 
571  CALL dlarrk( in, in, gl, gu, d(ibegin),
572  $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
573  IF( iinfo.NE.0 ) THEN
574  info = -1
575  RETURN
576  ENDIF
577  isrght = min(gu, tmp + tmp1
578  $ + hndrd * eps * abs(tmp + tmp1))
579 * Improve the estimate of the spectral diameter
580  spdiam = isrght - isleft
581  ELSE
582 * Case of bisection
583 * Find approximations to the wanted extremal eigenvalues
584  isleft = max(gl, w(wbegin) - werr(wbegin)
585  $ - hndrd * eps*abs(w(wbegin)- werr(wbegin) ))
586  isrght = min(gu,w(wend) + werr(wend)
587  $ + hndrd * eps * abs(w(wend)+ werr(wend)))
588  ENDIF
589 
590 
591 * Decide whether the base representation for the current block
592 * L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I
593 * should be on the left or the right end of the current block.
594 * The strategy is to shift to the end which is "more populated"
595 * Furthermore, decide whether to use DQDS for the computation of
596 * the eigenvalue approximations at the end of DLARRE or bisection.
597 * dqds is chosen if all eigenvalues are desired or the number of
598 * eigenvalues to be computed is large compared to the blocksize.
599  IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
600 * If all the eigenvalues have to be computed, we use dqd
601  usedqd = .true.
602 * INDL is the local index of the first eigenvalue to compute
603  indl = 1
604  indu = in
605 * MB = number of eigenvalues to compute
606  mb = in
607  wend = wbegin + mb - 1
608 * Define 1/4 and 3/4 points of the spectrum
609  s1 = isleft + fourth * spdiam
610  s2 = isrght - fourth * spdiam
611  ELSE
612 * DLARRD has computed IBLOCK and INDEXW for each eigenvalue
613 * approximation.
614 * choose sigma
615  IF( usedqd ) THEN
616  s1 = isleft + fourth * spdiam
617  s2 = isrght - fourth * spdiam
618  ELSE
619  tmp = min(isrght,vu) - max(isleft,vl)
620  s1 = max(isleft,vl) + fourth * tmp
621  s2 = min(isrght,vu) - fourth * tmp
622  ENDIF
623  ENDIF
624 
625 * Compute the negcount at the 1/4 and 3/4 points
626  IF(mb.GT.1) THEN
627  CALL dlarrc( 'T', in, s1, s2, d(ibegin),
628  $ e(ibegin), pivmin, cnt, cnt1, cnt2, iinfo)
629  ENDIF
630 
631  IF(mb.EQ.1) THEN
632  sigma = gl
633  sgndef = one
634  ELSEIF( cnt1 - indl .GE. indu - cnt2 ) THEN
635  IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
636  sigma = max(isleft,gl)
637  ELSEIF( usedqd ) THEN
638 * use Gerschgorin bound as shift to get pos def matrix
639 * for dqds
640  sigma = isleft
641  ELSE
642 * use approximation of the first desired eigenvalue of the
643 * block as shift
644  sigma = max(isleft,vl)
645  ENDIF
646  sgndef = one
647  ELSE
648  IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
649  sigma = min(isrght,gu)
650  ELSEIF( usedqd ) THEN
651 * use Gerschgorin bound as shift to get neg def matrix
652 * for dqds
653  sigma = isrght
654  ELSE
655 * use approximation of the first desired eigenvalue of the
656 * block as shift
657  sigma = min(isrght,vu)
658  ENDIF
659  sgndef = -one
660  ENDIF
661 
662 
663 * An initial SIGMA has been chosen that will be used for computing
664 * T - SIGMA I = L D L^T
665 * Define the increment TAU of the shift in case the initial shift
666 * needs to be refined to obtain a factorization with not too much
667 * element growth.
668  IF( usedqd ) THEN
669 * The initial SIGMA was to the outer end of the spectrum
670 * the matrix is definite and we need not retreat.
671  tau = spdiam*eps*n + two*pivmin
672  tau = max( tau,two*eps*abs(sigma) )
673  ELSE
674  IF(mb.GT.1) THEN
675  clwdth = w(wend) + werr(wend) - w(wbegin) - werr(wbegin)
676  avgap = abs(clwdth / dble(wend-wbegin))
677  IF( sgndef.EQ.one ) THEN
678  tau = half*max(wgap(wbegin),avgap)
679  tau = max(tau,werr(wbegin))
680  ELSE
681  tau = half*max(wgap(wend-1),avgap)
682  tau = max(tau,werr(wend))
683  ENDIF
684  ELSE
685  tau = werr(wbegin)
686  ENDIF
687  ENDIF
688 *
689  DO 80 idum = 1, maxtry
690 * Compute L D L^T factorization of tridiagonal matrix T - sigma I.
691 * Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
692 * pivots in WORK(2*IN+1:3*IN)
693  dpivot = d( ibegin ) - sigma
694  work( 1 ) = dpivot
695  dmax = abs( work(1) )
696  j = ibegin
697  DO 70 i = 1, in - 1
698  work( 2*in+i ) = one / work( i )
699  tmp = e( j )*work( 2*in+i )
700  work( in+i ) = tmp
701  dpivot = ( d( j+1 )-sigma ) - tmp*e( j )
702  work( i+1 ) = dpivot
703  dmax = max( dmax, abs(dpivot) )
704  j = j + 1
705  70 CONTINUE
706 * check for element growth
707  IF( dmax .GT. maxgrowth*spdiam ) THEN
708  norep = .true.
709  ELSE
710  norep = .false.
711  ENDIF
712  IF( usedqd .AND. .NOT.norep ) THEN
713 * Ensure the definiteness of the representation
714 * All entries of D (of L D L^T) must have the same sign
715  DO 71 i = 1, in
716  tmp = sgndef*work( i )
717  IF( tmp.LT.zero ) norep = .true.
718  71 CONTINUE
719  ENDIF
720  IF(norep) THEN
721 * Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin
722 * shift which makes the matrix definite. So we should end up
723 * here really only in the case of IRANGE = VALRNG or INDRNG.
724  IF( idum.EQ.maxtry-1 ) THEN
725  IF( sgndef.EQ.one ) THEN
726 * The fudged Gerschgorin shift should succeed
727  sigma =
728  $ gl - fudge*spdiam*eps*n - fudge*two*pivmin
729  ELSE
730  sigma =
731  $ gu + fudge*spdiam*eps*n + fudge*two*pivmin
732  END IF
733  ELSE
734  sigma = sigma - sgndef * tau
735  tau = two * tau
736  END IF
737  ELSE
738 * an initial RRR is found
739  GO TO 83
740  END IF
741  80 CONTINUE
742 * if the program reaches this point, no base representation could be
743 * found in MAXTRY iterations.
744  info = 2
745  RETURN
746 
747  83 CONTINUE
748 * At this point, we have found an initial base representation
749 * T - SIGMA I = L D L^T with not too much element growth.
750 * Store the shift.
751  e( iend ) = sigma
752 * Store D and L.
753  CALL dcopy( in, work, 1, d( ibegin ), 1 )
754  CALL dcopy( in-1, work( in+1 ), 1, e( ibegin ), 1 )
755 
756 
757  IF(mb.GT.1 ) THEN
758 *
759 * Perturb each entry of the base representation by a small
760 * (but random) relative amount to overcome difficulties with
761 * glued matrices.
762 *
763  DO 122 i = 1, 4
764  iseed( i ) = 1
765  122 CONTINUE
766 
767  CALL dlarnv(2, iseed, 2*in-1, work(1))
768  DO 125 i = 1,in-1
769  d(ibegin+i-1) = d(ibegin+i-1)*(one+eps*pert*work(i))
770  e(ibegin+i-1) = e(ibegin+i-1)*(one+eps*pert*work(in+i))
771  125 CONTINUE
772  d(iend) = d(iend)*(one+eps*four*work(in))
773 *
774  ENDIF
775 *
776 * Don't update the Gerschgorin intervals because keeping track
777 * of the updates would be too much work in DLARRV.
778 * We update W instead and use it to locate the proper Gerschgorin
779 * intervals.
780 
781 * Compute the required eigenvalues of L D L' by bisection or dqds
782  IF ( .NOT.usedqd ) THEN
783 * If DLARRD has been used, shift the eigenvalue approximations
784 * according to their representation. This is necessary for
785 * a uniform DLARRV since dqds computes eigenvalues of the
786 * shifted representation. In DLARRV, W will always hold the
787 * UNshifted eigenvalue approximation.
788  DO 134 j=wbegin,wend
789  w(j) = w(j) - sigma
790  werr(j) = werr(j) + abs(w(j)) * eps
791  134 CONTINUE
792 * call DLARRB to reduce eigenvalue error of the approximations
793 * from DLARRD
794  DO 135 i = ibegin, iend-1
795  work( i ) = d( i ) * e( i )**2
796  135 CONTINUE
797 * use bisection to find EV from INDL to INDU
798  CALL dlarrb(in, d(ibegin), work(ibegin),
799  $ indl, indu, rtol1, rtol2, indl-1,
800  $ w(wbegin), wgap(wbegin), werr(wbegin),
801  $ work( 2*n+1 ), iwork, pivmin, spdiam,
802  $ in, iinfo )
803  IF( iinfo .NE. 0 ) THEN
804  info = -4
805  RETURN
806  END IF
807 * DLARRB computes all gaps correctly except for the last one
808 * Record distance to VU/GU
809  wgap( wend ) = max( zero,
810  $ ( vu-sigma ) - ( w( wend ) + werr( wend ) ) )
811  DO 138 i = indl, indu
812  m = m + 1
813  iblock(m) = jblk
814  indexw(m) = i
815  138 CONTINUE
816  ELSE
817 * Call dqds to get all eigs (and then possibly delete unwanted
818 * eigenvalues).
819 * Note that dqds finds the eigenvalues of the L D L^T representation
820 * of T to high relative accuracy. High relative accuracy
821 * might be lost when the shift of the RRR is subtracted to obtain
822 * the eigenvalues of T. However, T is not guaranteed to define its
823 * eigenvalues to high relative accuracy anyway.
824 * Set RTOL to the order of the tolerance used in DLASQ2
825 * This is an ESTIMATED error, the worst case bound is 4*N*EPS
826 * which is usually too large and requires unnecessary work to be
827 * done by bisection when computing the eigenvectors
828  rtol = log(dble(in)) * four * eps
829  j = ibegin
830  DO 140 i = 1, in - 1
831  work( 2*i-1 ) = abs( d( j ) )
832  work( 2*i ) = e( j )*e( j )*work( 2*i-1 )
833  j = j + 1
834  140 CONTINUE
835  work( 2*in-1 ) = abs( d( iend ) )
836  work( 2*in ) = zero
837  CALL dlasq2( in, work, iinfo )
838  IF( iinfo .NE. 0 ) THEN
839 * If IINFO = -5 then an index is part of a tight cluster
840 * and should be changed. The index is in IWORK(1) and the
841 * gap is in WORK(N+1)
842  info = -5
843  RETURN
844  ELSE
845 * Test that all eigenvalues are positive as expected
846  DO 149 i = 1, in
847  IF( work( i ).LT.zero ) THEN
848  info = -6
849  RETURN
850  ENDIF
851  149 CONTINUE
852  END IF
853  IF( sgndef.GT.zero ) THEN
854  DO 150 i = indl, indu
855  m = m + 1
856  w( m ) = work( in-i+1 )
857  iblock( m ) = jblk
858  indexw( m ) = i
859  150 CONTINUE
860  ELSE
861  DO 160 i = indl, indu
862  m = m + 1
863  w( m ) = -work( i )
864  iblock( m ) = jblk
865  indexw( m ) = i
866  160 CONTINUE
867  END IF
868 
869  DO 165 i = m - mb + 1, m
870 * the value of RTOL below should be the tolerance in DLASQ2
871  werr( i ) = rtol * abs( w(i) )
872  165 CONTINUE
873  DO 166 i = m - mb + 1, m - 1
874 * compute the right gap between the intervals
875  wgap( i ) = max( zero,
876  $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
877  166 CONTINUE
878  wgap( m ) = max( zero,
879  $ ( vu-sigma ) - ( w( m ) + werr( m ) ) )
880  END IF
881 * proceed with next block
882  ibegin = iend + 1
883  wbegin = wend + 1
884  170 CONTINUE
885 *
886 
887  RETURN
888 *
889 * end of DLARRE
890 *
891  END
subroutine dlarrc(JOBT, N, VL, VU, D, E, PIVMIN, EIGCNT, LCNT, RCNT, INFO)
DLARRC computes the number of eigenvalues of the symmetric tridiagonal matrix.
Definition: dlarrc.f:138
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
subroutine dlarrb(N, D, LLD, IFIRST, ILAST, RTOL1, RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK, PIVMIN, SPDIAM, TWIST, INFO)
DLARRB provides limited bisection to locate eigenvalues for more accuracy.
Definition: dlarrb.f:198
subroutine dlarra(N, D, E, E2, SPLTOL, TNRM, NSPLIT, ISPLIT, INFO)
DLARRA computes the splitting points with the specified threshold.
Definition: dlarra.f:138
subroutine dlarrk(N, IW, GL, GU, D, E2, PIVMIN, RELTOL, W, WERR, INFO)
DLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy.
Definition: dlarrk.f:147
subroutine dlasq2(N, Z, INFO)
DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated ...
Definition: dlasq2.f:114
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: dlarnv.f:99
subroutine dlarre(RANGE, N, VL, VU, IL, IU, D, E, E2, RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M, W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN, WORK, IWORK, INFO)
DLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduce...
Definition: dlarre.f:299
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53