LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
slarre.f
Go to the documentation of this file.
1 *> \brief \b SLARRE 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 SLARRE + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarre.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarre.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarre.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLARRE( 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 * REAL PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
30 * ..
31 * .. Array Arguments ..
32 * INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ),
33 * $ INDEXW( * )
34 * REAL 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, SLARRE 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 *> SSTEMR 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 SLASQ2) to
54 *> conpute all and then discard any unwanted one.
55 *> As an added benefit, SLARRE 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 REAL
81 *> \endverbatim
82 *>
83 *> \param[in,out] VU
84 *> \verbatim
85 *> VU is REAL
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', SLARRE 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 REAL 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 REAL 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 REAL 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 REAL
138 *> \endverbatim
139 *>
140 *> \param[in] RTOL2
141 *> \verbatim
142 *> RTOL2 is REAL
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 REAL
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 REAL 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 ( SLARRE may use the
183 *> remaining N-M elements as workspace).
184 *> \endverbatim
185 *>
186 *> \param[out] WERR
187 *> \verbatim
188 *> WERR is REAL 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 REAL 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 REAL 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 REAL
228 *> The minimum pivot in the Sturm sequence for T.
229 *> \endverbatim
230 *>
231 *> \param[out] WORK
232 *> \verbatim
233 *> WORK is REAL 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 SLARRE.
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 SLARRD.
253 *> = 2: No base representation could be found in MAXTRY iterations.
254 *> Increasing MAXTRY and recompilation might be a remedy.
255 *> =-3: Problem in SLARRB when computing the refined root
256 *> representation for SLASQ2.
257 *> =-4: Problem in SLARRB when preforming bisection on the
258 *> desired part of the spectrum.
259 *> =-5: Problem in SLASQ2.
260 *> =-6: Problem in SLASQ2.
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 slarre( 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  REAL pivmin, rtol1, rtol2, spltol, vl, vu
309 * ..
310 * .. Array Arguments ..
311  INTEGER iblock( * ), isplit( * ), iwork( * ),
312  $ indexw( * )
313  REAL d( * ), e( * ), e2( * ), gers( * ),
314  $ w( * ),werr( * ), wgap( * ), work( * )
315 * ..
316 *
317 * =====================================================================
318 *
319 * .. Parameters ..
320  REAL fac, four, fourth, fudge, half, hndrd,
321  $ maxgrowth, one, pert, two, zero
322  parameter( zero = 0.0e0, one = 1.0e0,
323  $ two = 2.0e0, four=4.0e0,
324  $ hndrd = 100.0e0,
325  $ pert = 4.0e0,
326  $ half = one/two, fourth = one/four, fac= half,
327  $ maxgrowth = 64.0e0, fudge = 2.0e0 )
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  REAL 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  REAL slamch
350  EXTERNAL slamch, lsame
351 
352 * ..
353 * .. External Subroutines ..
354  EXTERNAL scopy, slarnv, slarra, slarrb, slarrc, slarrd,
355  $ slasq2
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 = slamch( 'S' )
381  eps = slamch( 'P' )
382 
383 * Set parameters
384  rtl = hndrd*eps
385 * If one were ever to ask for less initial precision in BSRTOL,
386 * one should keep in mind that for the subset case, the extremal
387 * eigenvalues must be at least as accurate as the current setting
388 * (eigenvalues in the middle need not as much accuracy)
389  bsrtol = sqrt(eps)*(0.5e-3)
390 
391 * Treat case of 1x1 matrix for quick return
392  IF( n.EQ.1 ) THEN
393  IF( (irange.EQ.allrng).OR.
394  $ ((irange.EQ.valrng).AND.(d(1).GT.vl).AND.(d(1).LE.vu)).OR.
395  $ ((irange.EQ.indrng).AND.(il.EQ.1).AND.(iu.EQ.1)) ) THEN
396  m = 1
397  w(1) = d(1)
398 * The computation error of the eigenvalue is zero
399  werr(1) = zero
400  wgap(1) = zero
401  iblock( 1 ) = 1
402  indexw( 1 ) = 1
403  gers(1) = d( 1 )
404  gers(2) = d( 1 )
405  ENDIF
406 * store the shift for the initial RRR, which is zero in this case
407  e(1) = zero
408  return
409  END IF
410 
411 * General case: tridiagonal matrix of order > 1
412 *
413 * Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter.
414 * Compute maximum off-diagonal entry and pivmin.
415  gl = d(1)
416  gu = d(1)
417  eold = zero
418  emax = zero
419  e(n) = zero
420  DO 5 i = 1,n
421  werr(i) = zero
422  wgap(i) = zero
423  eabs = abs( e(i) )
424  IF( eabs .GE. emax ) THEN
425  emax = eabs
426  END IF
427  tmp1 = eabs + eold
428  gers( 2*i-1) = d(i) - tmp1
429  gl = min( gl, gers( 2*i - 1))
430  gers( 2*i ) = d(i) + tmp1
431  gu = max( gu, gers(2*i) )
432  eold = eabs
433  5 continue
434 * The minimum pivot allowed in the Sturm sequence for T
435  pivmin = safmin * max( one, emax**2 )
436 * Compute spectral diameter. The Gerschgorin bounds give an
437 * estimate that is wrong by at most a factor of SQRT(2)
438  spdiam = gu - gl
439 
440 * Compute splitting points
441  CALL slarra( n, d, e, e2, spltol, spdiam,
442  $ nsplit, isplit, iinfo )
443 
444 * Can force use of bisection instead of faster DQDS.
445 * Option left in the code for future multisection work.
446  forceb = .false.
447 
448 * Initialize USEDQD, DQDS should be used for ALLRNG unless someone
449 * explicitly wants bisection.
450  usedqd = (( irange.EQ.allrng ) .AND. (.NOT.forceb))
451 
452  IF( (irange.EQ.allrng) .AND. (.NOT. forceb) ) THEN
453 * Set interval [VL,VU] that contains all eigenvalues
454  vl = gl
455  vu = gu
456  ELSE
457 * We call SLARRD to find crude approximations to the eigenvalues
458 * in the desired range. In case IRANGE = INDRNG, we also obtain the
459 * interval (VL,VU] that contains all the wanted eigenvalues.
460 * An interval [LEFT,RIGHT] has converged if
461 * RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT))
462 * SLARRD needs a WORK of size 4*N, IWORK of size 3*N
463  CALL slarrd( range, 'B', n, vl, vu, il, iu, gers,
464  $ bsrtol, d, e, e2, pivmin, nsplit, isplit,
465  $ mm, w, werr, vl, vu, iblock, indexw,
466  $ work, iwork, iinfo )
467  IF( iinfo.NE.0 ) THEN
468  info = -1
469  return
470  ENDIF
471 * Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0
472  DO 14 i = mm+1,n
473  w( i ) = zero
474  werr( i ) = zero
475  iblock( i ) = 0
476  indexw( i ) = 0
477  14 continue
478  END IF
479 
480 
481 ***
482 * Loop over unreduced blocks
483  ibegin = 1
484  wbegin = 1
485  DO 170 jblk = 1, nsplit
486  iend = isplit( jblk )
487  in = iend - ibegin + 1
488 
489 * 1 X 1 block
490  IF( in.EQ.1 ) THEN
491  IF( (irange.EQ.allrng).OR.( (irange.EQ.valrng).AND.
492  $ ( d( ibegin ).GT.vl ).AND.( d( ibegin ).LE.vu ) )
493  $ .OR. ( (irange.EQ.indrng).AND.(iblock(wbegin).EQ.jblk))
494  $ ) THEN
495  m = m + 1
496  w( m ) = d( ibegin )
497  werr(m) = zero
498 * The gap for a single block doesn't matter for the later
499 * algorithm and is assigned an arbitrary large value
500  wgap(m) = zero
501  iblock( m ) = jblk
502  indexw( m ) = 1
503  wbegin = wbegin + 1
504  ENDIF
505 * E( IEND ) holds the shift for the initial RRR
506  e( iend ) = zero
507  ibegin = iend + 1
508  go to 170
509  END IF
510 *
511 * Blocks of size larger than 1x1
512 *
513 * E( IEND ) will hold the shift for the initial RRR, for now set it =0
514  e( iend ) = zero
515 *
516 * Find local outer bounds GL,GU for the block
517  gl = d(ibegin)
518  gu = d(ibegin)
519  DO 15 i = ibegin , iend
520  gl = min( gers( 2*i-1 ), gl )
521  gu = max( gers( 2*i ), gu )
522  15 continue
523  spdiam = gu - gl
524 
525  IF(.NOT. ((irange.EQ.allrng).AND.(.NOT.forceb)) ) THEN
526 * Count the number of eigenvalues in the current block.
527  mb = 0
528  DO 20 i = wbegin,mm
529  IF( iblock(i).EQ.jblk ) THEN
530  mb = mb+1
531  ELSE
532  goto 21
533  ENDIF
534  20 continue
535  21 continue
536 
537  IF( mb.EQ.0) THEN
538 * No eigenvalue in the current block lies in the desired range
539 * E( IEND ) holds the shift for the initial RRR
540  e( iend ) = zero
541  ibegin = iend + 1
542  go to 170
543  ELSE
544 
545 * Decide whether dqds or bisection is more efficient
546  usedqd = ( (mb .GT. fac*in) .AND. (.NOT.forceb) )
547  wend = wbegin + mb - 1
548 * Calculate gaps for the current block
549 * In later stages, when representations for individual
550 * eigenvalues are different, we use SIGMA = E( IEND ).
551  sigma = zero
552  DO 30 i = wbegin, wend - 1
553  wgap( i ) = max( zero,
554  $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
555  30 continue
556  wgap( wend ) = max( zero,
557  $ vu - sigma - (w( wend )+werr( wend )))
558 * Find local index of the first and last desired evalue.
559  indl = indexw(wbegin)
560  indu = indexw( wend )
561  ENDIF
562  ENDIF
563  IF(( (irange.EQ.allrng) .AND. (.NOT. forceb) ).OR.usedqd) THEN
564 * Case of DQDS
565 * Find approximations to the extremal eigenvalues of the block
566  CALL slarrk( in, 1, gl, gu, d(ibegin),
567  $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
568  IF( iinfo.NE.0 ) THEN
569  info = -1
570  return
571  ENDIF
572  isleft = max(gl, tmp - tmp1
573  $ - hndrd * eps* abs(tmp - tmp1))
574 
575  CALL slarrk( in, in, gl, gu, d(ibegin),
576  $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
577  IF( iinfo.NE.0 ) THEN
578  info = -1
579  return
580  ENDIF
581  isrght = min(gu, tmp + tmp1
582  $ + hndrd * eps * abs(tmp + tmp1))
583 * Improve the estimate of the spectral diameter
584  spdiam = isrght - isleft
585  ELSE
586 * Case of bisection
587 * Find approximations to the wanted extremal eigenvalues
588  isleft = max(gl, w(wbegin) - werr(wbegin)
589  $ - hndrd * eps*abs(w(wbegin)- werr(wbegin) ))
590  isrght = min(gu,w(wend) + werr(wend)
591  $ + hndrd * eps * abs(w(wend)+ werr(wend)))
592  ENDIF
593 
594 
595 * Decide whether the base representation for the current block
596 * L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I
597 * should be on the left or the right end of the current block.
598 * The strategy is to shift to the end which is "more populated"
599 * Furthermore, decide whether to use DQDS for the computation of
600 * the eigenvalue approximations at the end of SLARRE or bisection.
601 * dqds is chosen if all eigenvalues are desired or the number of
602 * eigenvalues to be computed is large compared to the blocksize.
603  IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
604 * If all the eigenvalues have to be computed, we use dqd
605  usedqd = .true.
606 * INDL is the local index of the first eigenvalue to compute
607  indl = 1
608  indu = in
609 * MB = number of eigenvalues to compute
610  mb = in
611  wend = wbegin + mb - 1
612 * Define 1/4 and 3/4 points of the spectrum
613  s1 = isleft + fourth * spdiam
614  s2 = isrght - fourth * spdiam
615  ELSE
616 * SLARRD has computed IBLOCK and INDEXW for each eigenvalue
617 * approximation.
618 * choose sigma
619  IF( usedqd ) THEN
620  s1 = isleft + fourth * spdiam
621  s2 = isrght - fourth * spdiam
622  ELSE
623  tmp = min(isrght,vu) - max(isleft,vl)
624  s1 = max(isleft,vl) + fourth * tmp
625  s2 = min(isrght,vu) - fourth * tmp
626  ENDIF
627  ENDIF
628 
629 * Compute the negcount at the 1/4 and 3/4 points
630  IF(mb.GT.1) THEN
631  CALL slarrc( 'T', in, s1, s2, d(ibegin),
632  $ e(ibegin), pivmin, cnt, cnt1, cnt2, iinfo)
633  ENDIF
634 
635  IF(mb.EQ.1) THEN
636  sigma = gl
637  sgndef = one
638  elseif( cnt1 - indl .GE. indu - cnt2 ) THEN
639  IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
640  sigma = max(isleft,gl)
641  elseif( usedqd ) THEN
642 * use Gerschgorin bound as shift to get pos def matrix
643 * for dqds
644  sigma = isleft
645  ELSE
646 * use approximation of the first desired eigenvalue of the
647 * block as shift
648  sigma = max(isleft,vl)
649  ENDIF
650  sgndef = one
651  ELSE
652  IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
653  sigma = min(isrght,gu)
654  elseif( usedqd ) THEN
655 * use Gerschgorin bound as shift to get neg def matrix
656 * for dqds
657  sigma = isrght
658  ELSE
659 * use approximation of the first desired eigenvalue of the
660 * block as shift
661  sigma = min(isrght,vu)
662  ENDIF
663  sgndef = -one
664  ENDIF
665 
666 
667 * An initial SIGMA has been chosen that will be used for computing
668 * T - SIGMA I = L D L^T
669 * Define the increment TAU of the shift in case the initial shift
670 * needs to be refined to obtain a factorization with not too much
671 * element growth.
672  IF( usedqd ) THEN
673 * The initial SIGMA was to the outer end of the spectrum
674 * the matrix is definite and we need not retreat.
675  tau = spdiam*eps*n + two*pivmin
676  tau = max( tau,two*eps*abs(sigma) )
677  ELSE
678  IF(mb.GT.1) THEN
679  clwdth = w(wend) + werr(wend) - w(wbegin) - werr(wbegin)
680  avgap = abs(clwdth / REAL(wend-wbegin))
681  IF( sgndef.EQ.one ) THEN
682  tau = half*max(wgap(wbegin),avgap)
683  tau = max(tau,werr(wbegin))
684  ELSE
685  tau = half*max(wgap(wend-1),avgap)
686  tau = max(tau,werr(wend))
687  ENDIF
688  ELSE
689  tau = werr(wbegin)
690  ENDIF
691  ENDIF
692 *
693  DO 80 idum = 1, maxtry
694 * Compute L D L^T factorization of tridiagonal matrix T - sigma I.
695 * Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
696 * pivots in WORK(2*IN+1:3*IN)
697  dpivot = d( ibegin ) - sigma
698  work( 1 ) = dpivot
699  dmax = abs( work(1) )
700  j = ibegin
701  DO 70 i = 1, in - 1
702  work( 2*in+i ) = one / work( i )
703  tmp = e( j )*work( 2*in+i )
704  work( in+i ) = tmp
705  dpivot = ( d( j+1 )-sigma ) - tmp*e( j )
706  work( i+1 ) = dpivot
707  dmax = max( dmax, abs(dpivot) )
708  j = j + 1
709  70 continue
710 * check for element growth
711  IF( dmax .GT. maxgrowth*spdiam ) THEN
712  norep = .true.
713  ELSE
714  norep = .false.
715  ENDIF
716  IF( usedqd .AND. .NOT.norep ) THEN
717 * Ensure the definiteness of the representation
718 * All entries of D (of L D L^T) must have the same sign
719  DO 71 i = 1, in
720  tmp = sgndef*work( i )
721  IF( tmp.LT.zero ) norep = .true.
722  71 continue
723  ENDIF
724  IF(norep) THEN
725 * Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin
726 * shift which makes the matrix definite. So we should end up
727 * here really only in the case of IRANGE = VALRNG or INDRNG.
728  IF( idum.EQ.maxtry-1 ) THEN
729  IF( sgndef.EQ.one ) THEN
730 * The fudged Gerschgorin shift should succeed
731  sigma =
732  $ gl - fudge*spdiam*eps*n - fudge*two*pivmin
733  ELSE
734  sigma =
735  $ gu + fudge*spdiam*eps*n + fudge*two*pivmin
736  END IF
737  ELSE
738  sigma = sigma - sgndef * tau
739  tau = two * tau
740  END IF
741  ELSE
742 * an initial RRR is found
743  go to 83
744  END IF
745  80 continue
746 * if the program reaches this point, no base representation could be
747 * found in MAXTRY iterations.
748  info = 2
749  return
750 
751  83 continue
752 * At this point, we have found an initial base representation
753 * T - SIGMA I = L D L^T with not too much element growth.
754 * Store the shift.
755  e( iend ) = sigma
756 * Store D and L.
757  CALL scopy( in, work, 1, d( ibegin ), 1 )
758  CALL scopy( in-1, work( in+1 ), 1, e( ibegin ), 1 )
759 
760 
761  IF(mb.GT.1 ) THEN
762 *
763 * Perturb each entry of the base representation by a small
764 * (but random) relative amount to overcome difficulties with
765 * glued matrices.
766 *
767  DO 122 i = 1, 4
768  iseed( i ) = 1
769  122 continue
770 
771  CALL slarnv(2, iseed, 2*in-1, work(1))
772  DO 125 i = 1,in-1
773  d(ibegin+i-1) = d(ibegin+i-1)*(one+eps*pert*work(i))
774  e(ibegin+i-1) = e(ibegin+i-1)*(one+eps*pert*work(in+i))
775  125 continue
776  d(iend) = d(iend)*(one+eps*four*work(in))
777 *
778  ENDIF
779 *
780 * Don't update the Gerschgorin intervals because keeping track
781 * of the updates would be too much work in SLARRV.
782 * We update W instead and use it to locate the proper Gerschgorin
783 * intervals.
784 
785 * Compute the required eigenvalues of L D L' by bisection or dqds
786  IF ( .NOT.usedqd ) THEN
787 * If SLARRD has been used, shift the eigenvalue approximations
788 * according to their representation. This is necessary for
789 * a uniform SLARRV since dqds computes eigenvalues of the
790 * shifted representation. In SLARRV, W will always hold the
791 * UNshifted eigenvalue approximation.
792  DO 134 j=wbegin,wend
793  w(j) = w(j) - sigma
794  werr(j) = werr(j) + abs(w(j)) * eps
795  134 continue
796 * call SLARRB to reduce eigenvalue error of the approximations
797 * from SLARRD
798  DO 135 i = ibegin, iend-1
799  work( i ) = d( i ) * e( i )**2
800  135 continue
801 * use bisection to find EV from INDL to INDU
802  CALL slarrb(in, d(ibegin), work(ibegin),
803  $ indl, indu, rtol1, rtol2, indl-1,
804  $ w(wbegin), wgap(wbegin), werr(wbegin),
805  $ work( 2*n+1 ), iwork, pivmin, spdiam,
806  $ in, iinfo )
807  IF( iinfo .NE. 0 ) THEN
808  info = -4
809  return
810  END IF
811 * SLARRB computes all gaps correctly except for the last one
812 * Record distance to VU/GU
813  wgap( wend ) = max( zero,
814  $ ( vu-sigma ) - ( w( wend ) + werr( wend ) ) )
815  DO 138 i = indl, indu
816  m = m + 1
817  iblock(m) = jblk
818  indexw(m) = i
819  138 continue
820  ELSE
821 * Call dqds to get all eigs (and then possibly delete unwanted
822 * eigenvalues).
823 * Note that dqds finds the eigenvalues of the L D L^T representation
824 * of T to high relative accuracy. High relative accuracy
825 * might be lost when the shift of the RRR is subtracted to obtain
826 * the eigenvalues of T. However, T is not guaranteed to define its
827 * eigenvalues to high relative accuracy anyway.
828 * Set RTOL to the order of the tolerance used in SLASQ2
829 * This is an ESTIMATED error, the worst case bound is 4*N*EPS
830 * which is usually too large and requires unnecessary work to be
831 * done by bisection when computing the eigenvectors
832  rtol = log(REAL(in)) * four * eps
833  j = ibegin
834  DO 140 i = 1, in - 1
835  work( 2*i-1 ) = abs( d( j ) )
836  work( 2*i ) = e( j )*e( j )*work( 2*i-1 )
837  j = j + 1
838  140 continue
839  work( 2*in-1 ) = abs( d( iend ) )
840  work( 2*in ) = zero
841  CALL slasq2( in, work, iinfo )
842  IF( iinfo .NE. 0 ) THEN
843 * If IINFO = -5 then an index is part of a tight cluster
844 * and should be changed. The index is in IWORK(1) and the
845 * gap is in WORK(N+1)
846  info = -5
847  return
848  ELSE
849 * Test that all eigenvalues are positive as expected
850  DO 149 i = 1, in
851  IF( work( i ).LT.zero ) THEN
852  info = -6
853  return
854  ENDIF
855  149 continue
856  END IF
857  IF( sgndef.GT.zero ) THEN
858  DO 150 i = indl, indu
859  m = m + 1
860  w( m ) = work( in-i+1 )
861  iblock( m ) = jblk
862  indexw( m ) = i
863  150 continue
864  ELSE
865  DO 160 i = indl, indu
866  m = m + 1
867  w( m ) = -work( i )
868  iblock( m ) = jblk
869  indexw( m ) = i
870  160 continue
871  END IF
872 
873  DO 165 i = m - mb + 1, m
874 * the value of RTOL below should be the tolerance in SLASQ2
875  werr( i ) = rtol * abs( w(i) )
876  165 continue
877  DO 166 i = m - mb + 1, m - 1
878 * compute the right gap between the intervals
879  wgap( i ) = max( zero,
880  $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
881  166 continue
882  wgap( m ) = max( zero,
883  $ ( vu-sigma ) - ( w( m ) + werr( m ) ) )
884  END IF
885 * proceed with next block
886  ibegin = iend + 1
887  wbegin = wend + 1
888  170 continue
889 *
890 
891  return
892 *
893 * end of SLARRE
894 *
895  END