LAPACK  3.6.0
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 *> \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
subroutine slasq2(N, Z, INFO)
SLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated ...
Definition: slasq2.f:114
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:323
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
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 slarra(N, D, E, E2, SPLTOL, TNRM, NSPLIT, ISPLIT, INFO)
SLARRA computes the splitting points with the specified threshold.
Definition: slarra.f:138
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:138
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 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:299
subroutine slarnv(IDIST, ISEED, N, X)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: slarnv.f:99