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