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