LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
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