LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
dlarre.f
Go to the documentation of this file.
1 *> \brief \b DLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduced block Ti, finds base representations and eigenvalues.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DLARRE + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarre.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarre.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarre.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLARRE( RANGE, N, VL, VU, IL, IU, D, E, E2,
22 * RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
23 * W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
24 * WORK, IWORK, INFO )
25 *
26 * .. Scalar Arguments ..
27 * CHARACTER RANGE
28 * INTEGER IL, INFO, IU, M, N, NSPLIT
29 * DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
30 * ..
31 * .. Array Arguments ..
32 * INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ),
33 * $ INDEXW( * )
34 * DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ),
35 * $ W( * ),WERR( * ), WGAP( * ), WORK( * )
36 * ..
37 *
38 *
39 *> \par Purpose:
40 * =============
41 *>
42 *> \verbatim
43 *>
44 *> To find the desired eigenvalues of a given real symmetric
45 *> tridiagonal matrix T, DLARRE sets any "small" off-diagonal
46 *> elements to zero, and for each unreduced block T_i, it finds
47 *> (a) a suitable shift at one end of the block's spectrum,
48 *> (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
49 *> (c) eigenvalues of each L_i D_i L_i^T.
50 *> The representations and eigenvalues found are then used by
51 *> DSTEMR to compute the eigenvectors of T.
52 *> The accuracy varies depending on whether bisection is used to
53 *> find a few eigenvalues or the dqds algorithm (subroutine DLASQ2) to
54 *> conpute all and then discard any unwanted one.
55 *> As an added benefit, DLARRE also outputs the n
56 *> Gerschgorin intervals for the matrices L_i D_i L_i^T.
57 *> \endverbatim
58 *
59 * Arguments:
60 * ==========
61 *
62 *> \param[in] RANGE
63 *> \verbatim
64 *> RANGE is CHARACTER*1
65 *> = 'A': ("All") all eigenvalues will be found.
66 *> = 'V': ("Value") all eigenvalues in the half-open interval
67 *> (VL, VU] will be found.
68 *> = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
69 *> entire matrix) will be found.
70 *> \endverbatim
71 *>
72 *> \param[in] N
73 *> \verbatim
74 *> N is INTEGER
75 *> The order of the matrix. N > 0.
76 *> \endverbatim
77 *>
78 *> \param[in,out] VL
79 *> \verbatim
80 *> VL is DOUBLE PRECISION
81 *> 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', DLARRE computes bounds on the desired
85 *> part of the spectrum.
86 *> \endverbatim
87 *>
88 *> \param[in,out] VU
89 *> \verbatim
90 *> VU is DOUBLE PRECISION
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', DLARRE 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
146 *> \endverbatim
147 *>
148 *> \param[in] RTOL2
149 *> \verbatim
150 *> RTOL2 is DOUBLE PRECISION
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 DOUBLE PRECISION
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 DOUBLE PRECISION 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 ( DLARRE may use the
191 *> remaining N-M elements as workspace).
192 *> \endverbatim
193 *>
194 *> \param[out] WERR
195 *> \verbatim
196 *> WERR is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
236 *> The minimum pivot in the Sturm sequence for T.
237 *> \endverbatim
238 *>
239 *> \param[out] WORK
240 *> \verbatim
241 *> WORK is DOUBLE PRECISION 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 DLARRE.
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 DLARRD.
261 *> = 2: No base representation could be found in MAXTRY iterations.
262 *> Increasing MAXTRY and recompilation might be a remedy.
263 *> =-3: Problem in DLARRB when computing the refined root
264 *> representation for DLASQ2.
265 *> =-4: Problem in DLARRB when preforming bisection on the
266 *> desired part of the spectrum.
267 *> =-5: Problem in DLASQ2.
268 *> =-6: Problem in DLASQ2.
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 auxOTHERauxiliary
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 dlarre( 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.6.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  DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
317 * ..
318 * .. Array Arguments ..
319  INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ),
320  $ indexw( * )
321  DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ),
322  $ w( * ),werr( * ), wgap( * ), work( * )
323 * ..
324 *
325 * =====================================================================
326 *
327 * .. Parameters ..
328  DOUBLE PRECISION FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
329  $ maxgrowth, one, pert, two, zero
330  parameter ( zero = 0.0d0, one = 1.0d0,
331  $ two = 2.0d0, four=4.0d0,
332  $ hndrd = 100.0d0,
333  $ pert = 8.0d0,
334  $ half = one/two, fourth = one/four, fac= half,
335  $ maxgrowth = 64.0d0, fudge = 2.0d0 )
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  DOUBLE PRECISION 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  DOUBLE PRECISION DLAMCH
358  EXTERNAL dlamch, lsame
359 
360 * ..
361 * .. External Subroutines ..
362  EXTERNAL dcopy, dlarnv, dlarra, dlarrb, dlarrc, dlarrd,
363  $ dlasq2
364 * ..
365 * .. Intrinsic Functions ..
366  INTRINSIC abs, max, min
367 
368 * ..
369 * .. Executable Statements ..
370 *
371 
372  info = 0
373 
374 *
375 * Decode RANGE
376 *
377  IF( lsame( range, 'A' ) ) THEN
378  irange = allrng
379  ELSE IF( lsame( range, 'V' ) ) THEN
380  irange = valrng
381  ELSE IF( lsame( range, 'I' ) ) THEN
382  irange = indrng
383  END IF
384 
385  m = 0
386 
387 * Get machine constants
388  safmin = dlamch( 'S' )
389  eps = dlamch( 'P' )
390 
391 * Set parameters
392  rtl = sqrt(eps)
393  bsrtol = sqrt(eps)
394 
395 * Treat case of 1x1 matrix for quick return
396  IF( n.EQ.1 ) THEN
397  IF( (irange.EQ.allrng).OR.
398  $ ((irange.EQ.valrng).AND.(d(1).GT.vl).AND.(d(1).LE.vu)).OR.
399  $ ((irange.EQ.indrng).AND.(il.EQ.1).AND.(iu.EQ.1)) ) THEN
400  m = 1
401  w(1) = d(1)
402 * The computation error of the eigenvalue is zero
403  werr(1) = zero
404  wgap(1) = zero
405  iblock( 1 ) = 1
406  indexw( 1 ) = 1
407  gers(1) = d( 1 )
408  gers(2) = d( 1 )
409  ENDIF
410 * store the shift for the initial RRR, which is zero in this case
411  e(1) = zero
412  RETURN
413  END IF
414 
415 * General case: tridiagonal matrix of order > 1
416 *
417 * Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter.
418 * Compute maximum off-diagonal entry and pivmin.
419  gl = d(1)
420  gu = d(1)
421  eold = zero
422  emax = zero
423  e(n) = zero
424  DO 5 i = 1,n
425  werr(i) = zero
426  wgap(i) = zero
427  eabs = abs( e(i) )
428  IF( eabs .GE. emax ) THEN
429  emax = eabs
430  END IF
431  tmp1 = eabs + eold
432  gers( 2*i-1) = d(i) - tmp1
433  gl = min( gl, gers( 2*i - 1))
434  gers( 2*i ) = d(i) + tmp1
435  gu = max( gu, gers(2*i) )
436  eold = eabs
437  5 CONTINUE
438 * The minimum pivot allowed in the Sturm sequence for T
439  pivmin = safmin * max( one, emax**2 )
440 * Compute spectral diameter. The Gerschgorin bounds give an
441 * estimate that is wrong by at most a factor of SQRT(2)
442  spdiam = gu - gl
443 
444 * Compute splitting points
445  CALL dlarra( n, d, e, e2, spltol, spdiam,
446  $ nsplit, isplit, iinfo )
447 
448 * Can force use of bisection instead of faster DQDS.
449 * Option left in the code for future multisection work.
450  forceb = .false.
451 
452 * Initialize USEDQD, DQDS should be used for ALLRNG unless someone
453 * explicitly wants bisection.
454  usedqd = (( irange.EQ.allrng ) .AND. (.NOT.forceb))
455 
456  IF( (irange.EQ.allrng) .AND. (.NOT. forceb) ) THEN
457 * Set interval [VL,VU] that contains all eigenvalues
458  vl = gl
459  vu = gu
460  ELSE
461 * We call DLARRD to find crude approximations to the eigenvalues
462 * in the desired range. In case IRANGE = INDRNG, we also obtain the
463 * interval (VL,VU] that contains all the wanted eigenvalues.
464 * An interval [LEFT,RIGHT] has converged if
465 * RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT))
466 * DLARRD needs a WORK of size 4*N, IWORK of size 3*N
467  CALL dlarrd( range, 'B', n, vl, vu, il, iu, gers,
468  $ bsrtol, d, e, e2, pivmin, nsplit, isplit,
469  $ mm, w, werr, vl, vu, iblock, indexw,
470  $ work, iwork, iinfo )
471  IF( iinfo.NE.0 ) THEN
472  info = -1
473  RETURN
474  ENDIF
475 * Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0
476  DO 14 i = mm+1,n
477  w( i ) = zero
478  werr( i ) = zero
479  iblock( i ) = 0
480  indexw( i ) = 0
481  14 CONTINUE
482  END IF
483 
484 
485 ***
486 * Loop over unreduced blocks
487  ibegin = 1
488  wbegin = 1
489  DO 170 jblk = 1, nsplit
490  iend = isplit( jblk )
491  in = iend - ibegin + 1
492 
493 * 1 X 1 block
494  IF( in.EQ.1 ) THEN
495  IF( (irange.EQ.allrng).OR.( (irange.EQ.valrng).AND.
496  $ ( d( ibegin ).GT.vl ).AND.( d( ibegin ).LE.vu ) )
497  $ .OR. ( (irange.EQ.indrng).AND.(iblock(wbegin).EQ.jblk))
498  $ ) THEN
499  m = m + 1
500  w( m ) = d( ibegin )
501  werr(m) = zero
502 * The gap for a single block doesn't matter for the later
503 * algorithm and is assigned an arbitrary large value
504  wgap(m) = zero
505  iblock( m ) = jblk
506  indexw( m ) = 1
507  wbegin = wbegin + 1
508  ENDIF
509 * E( IEND ) holds the shift for the initial RRR
510  e( iend ) = zero
511  ibegin = iend + 1
512  GO TO 170
513  END IF
514 *
515 * Blocks of size larger than 1x1
516 *
517 * E( IEND ) will hold the shift for the initial RRR, for now set it =0
518  e( iend ) = zero
519 *
520 * Find local outer bounds GL,GU for the block
521  gl = d(ibegin)
522  gu = d(ibegin)
523  DO 15 i = ibegin , iend
524  gl = min( gers( 2*i-1 ), gl )
525  gu = max( gers( 2*i ), gu )
526  15 CONTINUE
527  spdiam = gu - gl
528 
529  IF(.NOT. ((irange.EQ.allrng).AND.(.NOT.forceb)) ) THEN
530 * Count the number of eigenvalues in the current block.
531  mb = 0
532  DO 20 i = wbegin,mm
533  IF( iblock(i).EQ.jblk ) THEN
534  mb = mb+1
535  ELSE
536  GOTO 21
537  ENDIF
538  20 CONTINUE
539  21 CONTINUE
540 
541  IF( mb.EQ.0) THEN
542 * No eigenvalue in the current block lies in the desired range
543 * E( IEND ) holds the shift for the initial RRR
544  e( iend ) = zero
545  ibegin = iend + 1
546  GO TO 170
547  ELSE
548 
549 * Decide whether dqds or bisection is more efficient
550  usedqd = ( (mb .GT. fac*in) .AND. (.NOT.forceb) )
551  wend = wbegin + mb - 1
552 * Calculate gaps for the current block
553 * In later stages, when representations for individual
554 * eigenvalues are different, we use SIGMA = E( IEND ).
555  sigma = zero
556  DO 30 i = wbegin, wend - 1
557  wgap( i ) = max( zero,
558  $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
559  30 CONTINUE
560  wgap( wend ) = max( zero,
561  $ vu - sigma - (w( wend )+werr( wend )))
562 * Find local index of the first and last desired evalue.
563  indl = indexw(wbegin)
564  indu = indexw( wend )
565  ENDIF
566  ENDIF
567  IF(( (irange.EQ.allrng) .AND. (.NOT. forceb) ).OR.usedqd) THEN
568 * Case of DQDS
569 * Find approximations to the extremal eigenvalues of the block
570  CALL dlarrk( in, 1, gl, gu, d(ibegin),
571  $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
572  IF( iinfo.NE.0 ) THEN
573  info = -1
574  RETURN
575  ENDIF
576  isleft = max(gl, tmp - tmp1
577  $ - hndrd * eps* abs(tmp - tmp1))
578 
579  CALL dlarrk( in, in, 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  isrght = min(gu, tmp + tmp1
586  $ + hndrd * eps * abs(tmp + tmp1))
587 * Improve the estimate of the spectral diameter
588  spdiam = isrght - isleft
589  ELSE
590 * Case of bisection
591 * Find approximations to the wanted extremal eigenvalues
592  isleft = max(gl, w(wbegin) - werr(wbegin)
593  $ - hndrd * eps*abs(w(wbegin)- werr(wbegin) ))
594  isrght = min(gu,w(wend) + werr(wend)
595  $ + hndrd * eps * abs(w(wend)+ werr(wend)))
596  ENDIF
597 
598 
599 * Decide whether the base representation for the current block
600 * L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I
601 * should be on the left or the right end of the current block.
602 * The strategy is to shift to the end which is "more populated"
603 * Furthermore, decide whether to use DQDS for the computation of
604 * the eigenvalue approximations at the end of DLARRE or bisection.
605 * dqds is chosen if all eigenvalues are desired or the number of
606 * eigenvalues to be computed is large compared to the blocksize.
607  IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
608 * If all the eigenvalues have to be computed, we use dqd
609  usedqd = .true.
610 * INDL is the local index of the first eigenvalue to compute
611  indl = 1
612  indu = in
613 * MB = number of eigenvalues to compute
614  mb = in
615  wend = wbegin + mb - 1
616 * Define 1/4 and 3/4 points of the spectrum
617  s1 = isleft + fourth * spdiam
618  s2 = isrght - fourth * spdiam
619  ELSE
620 * DLARRD has computed IBLOCK and INDEXW for each eigenvalue
621 * approximation.
622 * choose sigma
623  IF( usedqd ) THEN
624  s1 = isleft + fourth * spdiam
625  s2 = isrght - fourth * spdiam
626  ELSE
627  tmp = min(isrght,vu) - max(isleft,vl)
628  s1 = max(isleft,vl) + fourth * tmp
629  s2 = min(isrght,vu) - fourth * tmp
630  ENDIF
631  ENDIF
632 
633 * Compute the negcount at the 1/4 and 3/4 points
634  IF(mb.GT.1) THEN
635  CALL dlarrc( 'T', in, s1, s2, d(ibegin),
636  $ e(ibegin), pivmin, cnt, cnt1, cnt2, iinfo)
637  ENDIF
638 
639  IF(mb.EQ.1) THEN
640  sigma = gl
641  sgndef = one
642  ELSEIF( cnt1 - indl .GE. indu - cnt2 ) THEN
643  IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
644  sigma = max(isleft,gl)
645  ELSEIF( usedqd ) THEN
646 * use Gerschgorin bound as shift to get pos def matrix
647 * for dqds
648  sigma = isleft
649  ELSE
650 * use approximation of the first desired eigenvalue of the
651 * block as shift
652  sigma = max(isleft,vl)
653  ENDIF
654  sgndef = one
655  ELSE
656  IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
657  sigma = min(isrght,gu)
658  ELSEIF( usedqd ) THEN
659 * use Gerschgorin bound as shift to get neg def matrix
660 * for dqds
661  sigma = isrght
662  ELSE
663 * use approximation of the first desired eigenvalue of the
664 * block as shift
665  sigma = min(isrght,vu)
666  ENDIF
667  sgndef = -one
668  ENDIF
669 
670 
671 * An initial SIGMA has been chosen that will be used for computing
672 * T - SIGMA I = L D L^T
673 * Define the increment TAU of the shift in case the initial shift
674 * needs to be refined to obtain a factorization with not too much
675 * element growth.
676  IF( usedqd ) THEN
677 * The initial SIGMA was to the outer end of the spectrum
678 * the matrix is definite and we need not retreat.
679  tau = spdiam*eps*n + two*pivmin
680  tau = max( tau,two*eps*abs(sigma) )
681  ELSE
682  IF(mb.GT.1) THEN
683  clwdth = w(wend) + werr(wend) - w(wbegin) - werr(wbegin)
684  avgap = abs(clwdth / dble(wend-wbegin))
685  IF( sgndef.EQ.one ) THEN
686  tau = half*max(wgap(wbegin),avgap)
687  tau = max(tau,werr(wbegin))
688  ELSE
689  tau = half*max(wgap(wend-1),avgap)
690  tau = max(tau,werr(wend))
691  ENDIF
692  ELSE
693  tau = werr(wbegin)
694  ENDIF
695  ENDIF
696 *
697  DO 80 idum = 1, maxtry
698 * Compute L D L^T factorization of tridiagonal matrix T - sigma I.
699 * Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
700 * pivots in WORK(2*IN+1:3*IN)
701  dpivot = d( ibegin ) - sigma
702  work( 1 ) = dpivot
703  dmax = abs( work(1) )
704  j = ibegin
705  DO 70 i = 1, in - 1
706  work( 2*in+i ) = one / work( i )
707  tmp = e( j )*work( 2*in+i )
708  work( in+i ) = tmp
709  dpivot = ( d( j+1 )-sigma ) - tmp*e( j )
710  work( i+1 ) = dpivot
711  dmax = max( dmax, abs(dpivot) )
712  j = j + 1
713  70 CONTINUE
714 * check for element growth
715  IF( dmax .GT. maxgrowth*spdiam ) THEN
716  norep = .true.
717  ELSE
718  norep = .false.
719  ENDIF
720  IF( usedqd .AND. .NOT.norep ) THEN
721 * Ensure the definiteness of the representation
722 * All entries of D (of L D L^T) must have the same sign
723  DO 71 i = 1, in
724  tmp = sgndef*work( i )
725  IF( tmp.LT.zero ) norep = .true.
726  71 CONTINUE
727  ENDIF
728  IF(norep) THEN
729 * Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin
730 * shift which makes the matrix definite. So we should end up
731 * here really only in the case of IRANGE = VALRNG or INDRNG.
732  IF( idum.EQ.maxtry-1 ) THEN
733  IF( sgndef.EQ.one ) THEN
734 * The fudged Gerschgorin shift should succeed
735  sigma =
736  $ gl - fudge*spdiam*eps*n - fudge*two*pivmin
737  ELSE
738  sigma =
739  $ gu + fudge*spdiam*eps*n + fudge*two*pivmin
740  END IF
741  ELSE
742  sigma = sigma - sgndef * tau
743  tau = two * tau
744  END IF
745  ELSE
746 * an initial RRR is found
747  GO TO 83
748  END IF
749  80 CONTINUE
750 * if the program reaches this point, no base representation could be
751 * found in MAXTRY iterations.
752  info = 2
753  RETURN
754 
755  83 CONTINUE
756 * At this point, we have found an initial base representation
757 * T - SIGMA I = L D L^T with not too much element growth.
758 * Store the shift.
759  e( iend ) = sigma
760 * Store D and L.
761  CALL dcopy( in, work, 1, d( ibegin ), 1 )
762  CALL dcopy( in-1, work( in+1 ), 1, e( ibegin ), 1 )
763 
764 
765  IF(mb.GT.1 ) THEN
766 *
767 * Perturb each entry of the base representation by a small
768 * (but random) relative amount to overcome difficulties with
769 * glued matrices.
770 *
771  DO 122 i = 1, 4
772  iseed( i ) = 1
773  122 CONTINUE
774 
775  CALL dlarnv(2, iseed, 2*in-1, work(1))
776  DO 125 i = 1,in-1
777  d(ibegin+i-1) = d(ibegin+i-1)*(one+eps*pert*work(i))
778  e(ibegin+i-1) = e(ibegin+i-1)*(one+eps*pert*work(in+i))
779  125 CONTINUE
780  d(iend) = d(iend)*(one+eps*four*work(in))
781 *
782  ENDIF
783 *
784 * Don't update the Gerschgorin intervals because keeping track
785 * of the updates would be too much work in DLARRV.
786 * We update W instead and use it to locate the proper Gerschgorin
787 * intervals.
788 
789 * Compute the required eigenvalues of L D L' by bisection or dqds
790  IF ( .NOT.usedqd ) THEN
791 * If DLARRD has been used, shift the eigenvalue approximations
792 * according to their representation. This is necessary for
793 * a uniform DLARRV since dqds computes eigenvalues of the
794 * shifted representation. In DLARRV, W will always hold the
795 * UNshifted eigenvalue approximation.
796  DO 134 j=wbegin,wend
797  w(j) = w(j) - sigma
798  werr(j) = werr(j) + abs(w(j)) * eps
799  134 CONTINUE
800 * call DLARRB to reduce eigenvalue error of the approximations
801 * from DLARRD
802  DO 135 i = ibegin, iend-1
803  work( i ) = d( i ) * e( i )**2
804  135 CONTINUE
805 * use bisection to find EV from INDL to INDU
806  CALL dlarrb(in, d(ibegin), work(ibegin),
807  $ indl, indu, rtol1, rtol2, indl-1,
808  $ w(wbegin), wgap(wbegin), werr(wbegin),
809  $ work( 2*n+1 ), iwork, pivmin, spdiam,
810  $ in, iinfo )
811  IF( iinfo .NE. 0 ) THEN
812  info = -4
813  RETURN
814  END IF
815 * DLARRB computes all gaps correctly except for the last one
816 * Record distance to VU/GU
817  wgap( wend ) = max( zero,
818  $ ( vu-sigma ) - ( w( wend ) + werr( wend ) ) )
819  DO 138 i = indl, indu
820  m = m + 1
821  iblock(m) = jblk
822  indexw(m) = i
823  138 CONTINUE
824  ELSE
825 * Call dqds to get all eigs (and then possibly delete unwanted
826 * eigenvalues).
827 * Note that dqds finds the eigenvalues of the L D L^T representation
828 * of T to high relative accuracy. High relative accuracy
829 * might be lost when the shift of the RRR is subtracted to obtain
830 * the eigenvalues of T. However, T is not guaranteed to define its
831 * eigenvalues to high relative accuracy anyway.
832 * Set RTOL to the order of the tolerance used in DLASQ2
833 * This is an ESTIMATED error, the worst case bound is 4*N*EPS
834 * which is usually too large and requires unnecessary work to be
835 * done by bisection when computing the eigenvectors
836  rtol = log(dble(in)) * four * eps
837  j = ibegin
838  DO 140 i = 1, in - 1
839  work( 2*i-1 ) = abs( d( j ) )
840  work( 2*i ) = e( j )*e( j )*work( 2*i-1 )
841  j = j + 1
842  140 CONTINUE
843  work( 2*in-1 ) = abs( d( iend ) )
844  work( 2*in ) = zero
845  CALL dlasq2( in, work, iinfo )
846  IF( iinfo .NE. 0 ) THEN
847 * If IINFO = -5 then an index is part of a tight cluster
848 * and should be changed. The index is in IWORK(1) and the
849 * gap is in WORK(N+1)
850  info = -5
851  RETURN
852  ELSE
853 * Test that all eigenvalues are positive as expected
854  DO 149 i = 1, in
855  IF( work( i ).LT.zero ) THEN
856  info = -6
857  RETURN
858  ENDIF
859  149 CONTINUE
860  END IF
861  IF( sgndef.GT.zero ) THEN
862  DO 150 i = indl, indu
863  m = m + 1
864  w( m ) = work( in-i+1 )
865  iblock( m ) = jblk
866  indexw( m ) = i
867  150 CONTINUE
868  ELSE
869  DO 160 i = indl, indu
870  m = m + 1
871  w( m ) = -work( i )
872  iblock( m ) = jblk
873  indexw( m ) = i
874  160 CONTINUE
875  END IF
876 
877  DO 165 i = m - mb + 1, m
878 * the value of RTOL below should be the tolerance in DLASQ2
879  werr( i ) = rtol * abs( w(i) )
880  165 CONTINUE
881  DO 166 i = m - mb + 1, m - 1
882 * compute the right gap between the intervals
883  wgap( i ) = max( zero,
884  $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
885  166 CONTINUE
886  wgap( m ) = max( zero,
887  $ ( vu-sigma ) - ( w( m ) + werr( m ) ) )
888  END IF
889 * proceed with next block
890  ibegin = iend + 1
891  wbegin = wend + 1
892  170 CONTINUE
893 *
894 
895  RETURN
896 *
897 * end of DLARRE
898 *
899  END
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
subroutine dlarra(N, D, E, E2, SPLTOL, TNRM, NSPLIT, ISPLIT, INFO)
DLARRA computes the splitting points with the specified threshold.
Definition: dlarra.f:138
subroutine dlasq2(N, Z, INFO)
DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated ...
Definition: dlasq2.f:114
subroutine dlarrc(JOBT, N, VL, VU, D, E, PIVMIN, EIGCNT, LCNT, RCNT, INFO)
DLARRC computes the number of eigenvalues of the symmetric tridiagonal matrix.
Definition: dlarrc.f:139
subroutine dlarrb(N, D, LLD, IFIRST, ILAST, RTOL1, RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK, PIVMIN, SPDIAM, TWIST, INFO)
DLARRB provides limited bisection to locate eigenvalues for more accuracy.
Definition: dlarrb.f:198
subroutine dlarrk(N, IW, GL, GU, D, E2, PIVMIN, RELTOL, W, WERR, INFO)
DLARRK computes one eigenvalue of a symmetric tridiagonal matrix T to suitable accuracy.
Definition: dlarrk.f:147
subroutine dlarrd(RANGE, ORDER, N, VL, VU, IL, IU, GERS, RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, M, W, WERR, WL, WU, IBLOCK, INDEXW, WORK, IWORK, INFO)
DLARRD computes the eigenvalues of a symmetric tridiagonal matrix to suitable accuracy.
Definition: dlarrd.f:331
subroutine dlarre(RANGE, N, VL, VU, IL, IU, D, E, E2, RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M, W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN, WORK, IWORK, INFO)
DLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduce...
Definition: dlarre.f:307
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: dlarnv.f:99