LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> compute 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 < 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*> \ingroup larre
280*
281*> \par Further Details:
282* =====================
283*>
284*> \verbatim
285*>
286*> The base representations are required to suffer very little
287*> element growth and consequently define all their eigenvalues to
288*> high relative accuracy.
289*> \endverbatim
290*
291*> \par Contributors:
292* ==================
293*>
294*> Beresford Parlett, University of California, Berkeley, USA \n
295*> Jim Demmel, University of California, Berkeley, USA \n
296*> Inderjit Dhillon, University of Texas, Austin, USA \n
297*> Osni Marques, LBNL/NERSC, USA \n
298*> Christof Voemel, University of California, Berkeley, USA \n
299*>
300* =====================================================================
301 SUBROUTINE dlarre( RANGE, N, VL, VU, IL, IU, D, E, E2,
302 $ RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, M,
303 $ W, WERR, WGAP, IBLOCK, INDEXW, GERS, PIVMIN,
304 $ WORK, IWORK, INFO )
305*
306* -- LAPACK auxiliary routine --
307* -- LAPACK is a software package provided by Univ. of Tennessee, --
308* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
309*
310* .. Scalar Arguments ..
311 CHARACTER RANGE
312 INTEGER IL, INFO, IU, M, N, NSPLIT
313 DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU
314* ..
315* .. Array Arguments ..
316 INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ),
317 $ INDEXW( * )
318 DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ),
319 $ w( * ),werr( * ), wgap( * ), work( * )
320* ..
321*
322* =====================================================================
323*
324* .. Parameters ..
325 DOUBLE PRECISION FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD,
326 $ MAXGROWTH, ONE, PERT, TWO, ZERO
327 PARAMETER ( ZERO = 0.0d0, one = 1.0d0,
328 $ two = 2.0d0, four=4.0d0,
329 $ hndrd = 100.0d0,
330 $ pert = 8.0d0,
331 $ half = one/two, fourth = one/four, fac= half,
332 $ maxgrowth = 64.0d0, fudge = 2.0d0 )
333 INTEGER MAXTRY, ALLRNG, INDRNG, VALRNG
334 PARAMETER ( MAXTRY = 6, allrng = 1, indrng = 2,
335 $ valrng = 3 )
336* ..
337* .. Local Scalars ..
338 LOGICAL FORCEB, NOREP, USEDQD
339 INTEGER CNT, CNT1, CNT2, I, IBEGIN, IDUM, IEND, IINFO,
340 $ IN, INDL, INDU, IRANGE, J, JBLK, MB, MM,
341 $ wbegin, wend
342 DOUBLE PRECISION AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS,
343 $ EMAX, EOLD, EPS, GL, GU, ISLEFT, ISRGHT, RTL,
344 $ RTOL, S1, S2, SAFMIN, SGNDEF, SIGMA, SPDIAM,
345 $ tau, tmp, tmp1
346
347
348* ..
349* .. Local Arrays ..
350 INTEGER ISEED( 4 )
351* ..
352* .. External Functions ..
353 LOGICAL LSAME
354 DOUBLE PRECISION DLAMCH
355 EXTERNAL DLAMCH, LSAME
356
357* ..
358* .. External Subroutines ..
359 EXTERNAL dcopy, dlarnv, dlarra, dlarrb, dlarrc, dlarrd,
360 $ dlasq2, dlarrk
361* ..
362* .. Intrinsic Functions ..
363 INTRINSIC abs, max, min
364
365* ..
366* .. Executable Statements ..
367*
368
369 info = 0
370 nsplit = 0
371 m = 0
372*
373* Quick return if possible
374*
375 IF( n.LE.0 ) THEN
376 RETURN
377 END IF
378*
379* Decode RANGE
380*
381 IF( lsame( range, 'A' ) ) THEN
382 irange = allrng
383 ELSE IF( lsame( range, 'V' ) ) THEN
384 irange = valrng
385 ELSE IF( lsame( range, 'I' ) ) THEN
386 irange = indrng
387 END IF
388
389* Get machine constants
390 safmin = dlamch( 'S' )
391 eps = dlamch( 'P' )
392
393* Set parameters
394 rtl = sqrt(eps)
395 bsrtol = sqrt(eps)
396
397* Treat case of 1x1 matrix for quick return
398 IF( n.EQ.1 ) THEN
399 IF( (irange.EQ.allrng).OR.
400 $ ((irange.EQ.valrng).AND.(d(1).GT.vl).AND.(d(1).LE.vu)).OR.
401 $ ((irange.EQ.indrng).AND.(il.EQ.1).AND.(iu.EQ.1)) ) THEN
402 m = 1
403 w(1) = d(1)
404* The computation error of the eigenvalue is zero
405 werr(1) = zero
406 wgap(1) = zero
407 iblock( 1 ) = 1
408 indexw( 1 ) = 1
409 gers(1) = d( 1 )
410 gers(2) = d( 1 )
411 ENDIF
412* store the shift for the initial RRR, which is zero in this case
413 e(1) = zero
414 RETURN
415 END IF
416
417* General case: tridiagonal matrix of order > 1
418*
419* Init WERR, WGAP. Compute Gerschgorin intervals and spectral diameter.
420* Compute maximum off-diagonal entry and pivmin.
421 gl = d(1)
422 gu = d(1)
423 eold = zero
424 emax = zero
425 e(n) = zero
426 DO 5 i = 1,n
427 werr(i) = zero
428 wgap(i) = zero
429 eabs = abs( e(i) )
430 IF( eabs .GE. emax ) THEN
431 emax = eabs
432 END IF
433 tmp1 = eabs + eold
434 gers( 2*i-1) = d(i) - tmp1
435 gl = min( gl, gers( 2*i - 1))
436 gers( 2*i ) = d(i) + tmp1
437 gu = max( gu, gers(2*i) )
438 eold = eabs
439 5 CONTINUE
440* The minimum pivot allowed in the Sturm sequence for T
441 pivmin = safmin * max( one, emax**2 )
442* Compute spectral diameter. The Gerschgorin bounds give an
443* estimate that is wrong by at most a factor of SQRT(2)
444 spdiam = gu - gl
445
446* Compute splitting points
447 CALL dlarra( n, d, e, e2, spltol, spdiam,
448 $ nsplit, isplit, iinfo )
449
450* Can force use of bisection instead of faster DQDS.
451* Option left in the code for future multisection work.
452 forceb = .false.
453
454* Initialize USEDQD, DQDS should be used for ALLRNG unless someone
455* explicitly wants bisection.
456 usedqd = (( irange.EQ.allrng ) .AND. (.NOT.forceb))
457
458 IF( (irange.EQ.allrng) .AND. (.NOT. forceb) ) THEN
459* Set interval [VL,VU] that contains all eigenvalues
460 vl = gl
461 vu = gu
462 ELSE
463* We call DLARRD to find crude approximations to the eigenvalues
464* in the desired range. In case IRANGE = INDRNG, we also obtain the
465* interval (VL,VU] that contains all the wanted eigenvalues.
466* An interval [LEFT,RIGHT] has converged if
467* RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT))
468* DLARRD needs a WORK of size 4*N, IWORK of size 3*N
469 CALL dlarrd( range, 'B', n, vl, vu, il, iu, gers,
470 $ bsrtol, d, e, e2, pivmin, nsplit, isplit,
471 $ mm, w, werr, vl, vu, iblock, indexw,
472 $ work, iwork, iinfo )
473 IF( iinfo.NE.0 ) THEN
474 info = -1
475 RETURN
476 ENDIF
477* Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0
478 DO 14 i = mm+1,n
479 w( i ) = zero
480 werr( i ) = zero
481 iblock( i ) = 0
482 indexw( i ) = 0
483 14 CONTINUE
484 END IF
485
486
487***
488* Loop over unreduced blocks
489 ibegin = 1
490 wbegin = 1
491 DO 170 jblk = 1, nsplit
492 iend = isplit( jblk )
493 in = iend - ibegin + 1
494
495* 1 X 1 block
496 IF( in.EQ.1 ) THEN
497 IF( (irange.EQ.allrng).OR.( (irange.EQ.valrng).AND.
498 $ ( d( ibegin ).GT.vl ).AND.( d( ibegin ).LE.vu ) )
499 $ .OR. ( (irange.EQ.indrng).AND.(iblock(wbegin).EQ.jblk))
500 $ ) THEN
501 m = m + 1
502 w( m ) = d( ibegin )
503 werr(m) = zero
504* The gap for a single block doesn't matter for the later
505* algorithm and is assigned an arbitrary large value
506 wgap(m) = zero
507 iblock( m ) = jblk
508 indexw( m ) = 1
509 wbegin = wbegin + 1
510 ENDIF
511* E( IEND ) holds the shift for the initial RRR
512 e( iend ) = zero
513 ibegin = iend + 1
514 GO TO 170
515 END IF
516*
517* Blocks of size larger than 1x1
518*
519* E( IEND ) will hold the shift for the initial RRR, for now set it =0
520 e( iend ) = zero
521*
522* Find local outer bounds GL,GU for the block
523 gl = d(ibegin)
524 gu = d(ibegin)
525 DO 15 i = ibegin , iend
526 gl = min( gers( 2*i-1 ), gl )
527 gu = max( gers( 2*i ), gu )
528 15 CONTINUE
529 spdiam = gu - gl
530
531 IF(.NOT. ((irange.EQ.allrng).AND.(.NOT.forceb)) ) THEN
532* Count the number of eigenvalues in the current block.
533 mb = 0
534 DO 20 i = wbegin,mm
535 IF( iblock(i).EQ.jblk ) THEN
536 mb = mb+1
537 ELSE
538 GOTO 21
539 ENDIF
540 20 CONTINUE
541 21 CONTINUE
542
543 IF( mb.EQ.0) THEN
544* No eigenvalue in the current block lies in the desired range
545* E( IEND ) holds the shift for the initial RRR
546 e( iend ) = zero
547 ibegin = iend + 1
548 GO TO 170
549 ELSE
550
551* Decide whether dqds or bisection is more efficient
552 usedqd = ( (mb .GT. fac*in) .AND. (.NOT.forceb) )
553 wend = wbegin + mb - 1
554* Calculate gaps for the current block
555* In later stages, when representations for individual
556* eigenvalues are different, we use SIGMA = E( IEND ).
557 sigma = zero
558 DO 30 i = wbegin, wend - 1
559 wgap( i ) = max( zero,
560 $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
561 30 CONTINUE
562 wgap( wend ) = max( zero,
563 $ vu - sigma - (w( wend )+werr( wend )))
564* Find local index of the first and last desired evalue.
565 indl = indexw(wbegin)
566 indu = indexw( wend )
567 ENDIF
568 ENDIF
569 IF(( (irange.EQ.allrng) .AND. (.NOT. forceb) ).OR.usedqd) THEN
570* Case of DQDS
571* Find approximations to the extremal eigenvalues of the block
572 CALL dlarrk( in, 1, gl, gu, d(ibegin),
573 $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
574 IF( iinfo.NE.0 ) THEN
575 info = -1
576 RETURN
577 ENDIF
578 isleft = max(gl, tmp - tmp1
579 $ - hndrd * eps* abs(tmp - tmp1))
580
581 CALL dlarrk( in, in, gl, gu, d(ibegin),
582 $ e2(ibegin), pivmin, rtl, tmp, tmp1, iinfo )
583 IF( iinfo.NE.0 ) THEN
584 info = -1
585 RETURN
586 ENDIF
587 isrght = min(gu, tmp + tmp1
588 $ + hndrd * eps * abs(tmp + tmp1))
589* Improve the estimate of the spectral diameter
590 spdiam = isrght - isleft
591 ELSE
592* Case of bisection
593* Find approximations to the wanted extremal eigenvalues
594 isleft = max(gl, w(wbegin) - werr(wbegin)
595 $ - hndrd * eps*abs(w(wbegin)- werr(wbegin) ))
596 isrght = min(gu,w(wend) + werr(wend)
597 $ + hndrd * eps * abs(w(wend)+ werr(wend)))
598 ENDIF
599
600
601* Decide whether the base representation for the current block
602* L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I
603* should be on the left or the right end of the current block.
604* The strategy is to shift to the end which is "more populated"
605* Furthermore, decide whether to use DQDS for the computation of
606* the eigenvalue approximations at the end of DLARRE or bisection.
607* dqds is chosen if all eigenvalues are desired or the number of
608* eigenvalues to be computed is large compared to the blocksize.
609 IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
610* If all the eigenvalues have to be computed, we use dqd
611 usedqd = .true.
612* INDL is the local index of the first eigenvalue to compute
613 indl = 1
614 indu = in
615* MB = number of eigenvalues to compute
616 mb = in
617 wend = wbegin + mb - 1
618* Define 1/4 and 3/4 points of the spectrum
619 s1 = isleft + fourth * spdiam
620 s2 = isrght - fourth * spdiam
621 ELSE
622* DLARRD has computed IBLOCK and INDEXW for each eigenvalue
623* approximation.
624* choose sigma
625 IF( usedqd ) THEN
626 s1 = isleft + fourth * spdiam
627 s2 = isrght - fourth * spdiam
628 ELSE
629 tmp = min(isrght,vu) - max(isleft,vl)
630 s1 = max(isleft,vl) + fourth * tmp
631 s2 = min(isrght,vu) - fourth * tmp
632 ENDIF
633 ENDIF
634
635* Compute the negcount at the 1/4 and 3/4 points
636 IF(mb.GT.1) THEN
637 CALL dlarrc( 'T', in, s1, s2, d(ibegin),
638 $ e(ibegin), pivmin, cnt, cnt1, cnt2, iinfo)
639 ENDIF
640
641 IF(mb.EQ.1) THEN
642 sigma = gl
643 sgndef = one
644 ELSEIF( cnt1 - indl .GE. indu - cnt2 ) THEN
645 IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
646 sigma = max(isleft,gl)
647 ELSEIF( usedqd ) THEN
648* use Gerschgorin bound as shift to get pos def matrix
649* for dqds
650 sigma = isleft
651 ELSE
652* use approximation of the first desired eigenvalue of the
653* block as shift
654 sigma = max(isleft,vl)
655 ENDIF
656 sgndef = one
657 ELSE
658 IF( ( irange.EQ.allrng ) .AND. (.NOT.forceb) ) THEN
659 sigma = min(isrght,gu)
660 ELSEIF( usedqd ) THEN
661* use Gerschgorin bound as shift to get neg def matrix
662* for dqds
663 sigma = isrght
664 ELSE
665* use approximation of the first desired eigenvalue of the
666* block as shift
667 sigma = min(isrght,vu)
668 ENDIF
669 sgndef = -one
670 ENDIF
671
672
673* An initial SIGMA has been chosen that will be used for computing
674* T - SIGMA I = L D L^T
675* Define the increment TAU of the shift in case the initial shift
676* needs to be refined to obtain a factorization with not too much
677* element growth.
678 IF( usedqd ) THEN
679* The initial SIGMA was to the outer end of the spectrum
680* the matrix is definite and we need not retreat.
681 tau = spdiam*eps*n + two*pivmin
682 tau = max( tau,two*eps*abs(sigma) )
683 ELSE
684 IF(mb.GT.1) THEN
685 clwdth = w(wend) + werr(wend) - w(wbegin) - werr(wbegin)
686 avgap = abs(clwdth / dble(wend-wbegin))
687 IF( sgndef.EQ.one ) THEN
688 tau = half*max(wgap(wbegin),avgap)
689 tau = max(tau,werr(wbegin))
690 ELSE
691 tau = half*max(wgap(wend-1),avgap)
692 tau = max(tau,werr(wend))
693 ENDIF
694 ELSE
695 tau = werr(wbegin)
696 ENDIF
697 ENDIF
698*
699 DO 80 idum = 1, maxtry
700* Compute L D L^T factorization of tridiagonal matrix T - sigma I.
701* Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of
702* pivots in WORK(2*IN+1:3*IN)
703 dpivot = d( ibegin ) - sigma
704 work( 1 ) = dpivot
705 dmax = abs( work(1) )
706 j = ibegin
707 DO 70 i = 1, in - 1
708 work( 2*in+i ) = one / work( i )
709 tmp = e( j )*work( 2*in+i )
710 work( in+i ) = tmp
711 dpivot = ( d( j+1 )-sigma ) - tmp*e( j )
712 work( i+1 ) = dpivot
713 dmax = max( dmax, abs(dpivot) )
714 j = j + 1
715 70 CONTINUE
716* check for element growth
717 IF( dmax .GT. maxgrowth*spdiam ) THEN
718 norep = .true.
719 ELSE
720 norep = .false.
721 ENDIF
722 IF( usedqd .AND. .NOT.norep ) THEN
723* Ensure the definiteness of the representation
724* All entries of D (of L D L^T) must have the same sign
725 DO 71 i = 1, in
726 tmp = sgndef*work( i )
727 IF( tmp.LT.zero ) norep = .true.
728 71 CONTINUE
729 ENDIF
730 IF(norep) THEN
731* Note that in the case of IRANGE=ALLRNG, we use the Gerschgorin
732* shift which makes the matrix definite. So we should end up
733* here really only in the case of IRANGE = VALRNG or INDRNG.
734 IF( idum.EQ.maxtry-1 ) THEN
735 IF( sgndef.EQ.one ) THEN
736* The fudged Gerschgorin shift should succeed
737 sigma =
738 $ gl - fudge*spdiam*eps*n - fudge*two*pivmin
739 ELSE
740 sigma =
741 $ gu + fudge*spdiam*eps*n + fudge*two*pivmin
742 END IF
743 ELSE
744 sigma = sigma - sgndef * tau
745 tau = two * tau
746 END IF
747 ELSE
748* an initial RRR is found
749 GO TO 83
750 END IF
751 80 CONTINUE
752* if the program reaches this point, no base representation could be
753* found in MAXTRY iterations.
754 info = 2
755 RETURN
756
757 83 CONTINUE
758* At this point, we have found an initial base representation
759* T - SIGMA I = L D L^T with not too much element growth.
760* Store the shift.
761 e( iend ) = sigma
762* Store D and L.
763 CALL dcopy( in, work, 1, d( ibegin ), 1 )
764 CALL dcopy( in-1, work( in+1 ), 1, e( ibegin ), 1 )
765
766
767 IF(mb.GT.1 ) THEN
768*
769* Perturb each entry of the base representation by a small
770* (but random) relative amount to overcome difficulties with
771* glued matrices.
772*
773 DO 122 i = 1, 4
774 iseed( i ) = 1
775 122 CONTINUE
776
777 CALL dlarnv(2, iseed, 2*in-1, work(1))
778 DO 125 i = 1,in-1
779 d(ibegin+i-1) = d(ibegin+i-1)*(one+eps*pert*work(i))
780 e(ibegin+i-1) = e(ibegin+i-1)*(one+eps*pert*work(in+i))
781 125 CONTINUE
782 d(iend) = d(iend)*(one+eps*four*work(in))
783*
784 ENDIF
785*
786* Don't update the Gerschgorin intervals because keeping track
787* of the updates would be too much work in DLARRV.
788* We update W instead and use it to locate the proper Gerschgorin
789* intervals.
790
791* Compute the required eigenvalues of L D L' by bisection or dqds
792 IF ( .NOT.usedqd ) THEN
793* If DLARRD has been used, shift the eigenvalue approximations
794* according to their representation. This is necessary for
795* a uniform DLARRV since dqds computes eigenvalues of the
796* shifted representation. In DLARRV, W will always hold the
797* UNshifted eigenvalue approximation.
798 DO 134 j=wbegin,wend
799 w(j) = w(j) - sigma
800 werr(j) = werr(j) + abs(w(j)) * eps
801 134 CONTINUE
802* call DLARRB to reduce eigenvalue error of the approximations
803* from DLARRD
804 DO 135 i = ibegin, iend-1
805 work( i ) = d( i ) * e( i )**2
806 135 CONTINUE
807* use bisection to find EV from INDL to INDU
808 CALL dlarrb(in, d(ibegin), work(ibegin),
809 $ indl, indu, rtol1, rtol2, indl-1,
810 $ w(wbegin), wgap(wbegin), werr(wbegin),
811 $ work( 2*n+1 ), iwork, pivmin, spdiam,
812 $ in, iinfo )
813 IF( iinfo .NE. 0 ) THEN
814 info = -4
815 RETURN
816 END IF
817* DLARRB computes all gaps correctly except for the last one
818* Record distance to VU/GU
819 wgap( wend ) = max( zero,
820 $ ( vu-sigma ) - ( w( wend ) + werr( wend ) ) )
821 DO 138 i = indl, indu
822 m = m + 1
823 iblock(m) = jblk
824 indexw(m) = i
825 138 CONTINUE
826 ELSE
827* Call dqds to get all eigs (and then possibly delete unwanted
828* eigenvalues).
829* Note that dqds finds the eigenvalues of the L D L^T representation
830* of T to high relative accuracy. High relative accuracy
831* might be lost when the shift of the RRR is subtracted to obtain
832* the eigenvalues of T. However, T is not guaranteed to define its
833* eigenvalues to high relative accuracy anyway.
834* Set RTOL to the order of the tolerance used in DLASQ2
835* This is an ESTIMATED error, the worst case bound is 4*N*EPS
836* which is usually too large and requires unnecessary work to be
837* done by bisection when computing the eigenvectors
838 rtol = log(dble(in)) * four * eps
839 j = ibegin
840 DO 140 i = 1, in - 1
841 work( 2*i-1 ) = abs( d( j ) )
842 work( 2*i ) = e( j )*e( j )*work( 2*i-1 )
843 j = j + 1
844 140 CONTINUE
845 work( 2*in-1 ) = abs( d( iend ) )
846 work( 2*in ) = zero
847 CALL dlasq2( in, work, iinfo )
848 IF( iinfo .NE. 0 ) THEN
849* If IINFO = -5 then an index is part of a tight cluster
850* and should be changed. The index is in IWORK(1) and the
851* gap is in WORK(N+1)
852 info = -5
853 RETURN
854 ELSE
855* Test that all eigenvalues are positive as expected
856 DO 149 i = 1, in
857 IF( work( i ).LT.zero ) THEN
858 info = -6
859 RETURN
860 ENDIF
861 149 CONTINUE
862 END IF
863 IF( sgndef.GT.zero ) THEN
864 DO 150 i = indl, indu
865 m = m + 1
866 w( m ) = work( in-i+1 )
867 iblock( m ) = jblk
868 indexw( m ) = i
869 150 CONTINUE
870 ELSE
871 DO 160 i = indl, indu
872 m = m + 1
873 w( m ) = -work( i )
874 iblock( m ) = jblk
875 indexw( m ) = i
876 160 CONTINUE
877 END IF
878
879 DO 165 i = m - mb + 1, m
880* the value of RTOL below should be the tolerance in DLASQ2
881 werr( i ) = rtol * abs( w(i) )
882 165 CONTINUE
883 DO 166 i = m - mb + 1, m - 1
884* compute the right gap between the intervals
885 wgap( i ) = max( zero,
886 $ w(i+1)-werr(i+1) - (w(i)+werr(i)) )
887 166 CONTINUE
888 wgap( m ) = max( zero,
889 $ ( vu-sigma ) - ( w( m ) + werr( m ) ) )
890 END IF
891* proceed with next block
892 ibegin = iend + 1
893 wbegin = wend + 1
894 170 CONTINUE
895*
896
897 RETURN
898*
899* End of DLARRE
900*
901 END
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dlarnv(idist, iseed, n, x)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition dlarnv.f:97
subroutine dlarra(n, d, e, e2, spltol, tnrm, nsplit, isplit, info)
DLARRA computes the splitting points with the specified threshold.
Definition dlarra.f:136
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:196
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:137
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:329
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:305
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:145
subroutine dlasq2(n, z, info)
DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated ...
Definition dlasq2.f:112