LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlarre()

subroutine dlarre ( character  range,
integer  n,
double precision  vl,
double precision  vu,
integer  il,
integer  iu,
double precision, dimension( * )  d,
double precision, dimension( * )  e,
double precision, dimension( * )  e2,
double precision  rtol1,
double precision  rtol2,
double precision  spltol,
integer  nsplit,
integer, dimension( * )  isplit,
integer  m,
double precision, dimension( * )  w,
double precision, dimension( * )  werr,
double precision, dimension( * )  wgap,
integer, dimension( * )  iblock,
integer, dimension( * )  indexw,
double precision, dimension( * )  gers,
double precision  pivmin,
double precision, dimension( * )  work,
integer, dimension( * )  iwork,
integer  info 
)

DLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduced block Ti, finds base representations and eigenvalues.

Download DLARRE + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 To find the desired eigenvalues of a given real symmetric
 tridiagonal matrix T, DLARRE sets any "small" off-diagonal
 elements to zero, and for each unreduced block T_i, it finds
 (a) a suitable shift at one end of the block's spectrum,
 (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and
 (c) eigenvalues of each L_i D_i L_i^T.
 The representations and eigenvalues found are then used by
 DSTEMR to compute the eigenvectors of T.
 The accuracy varies depending on whether bisection is used to
 find a few eigenvalues or the dqds algorithm (subroutine DLASQ2) to
 compute all and then discard any unwanted one.
 As an added benefit, DLARRE also outputs the n
 Gerschgorin intervals for the matrices L_i D_i L_i^T.
Parameters
[in]RANGE
          RANGE is CHARACTER*1
          = 'A': ("All")   all eigenvalues will be found.
          = 'V': ("Value") all eigenvalues in the half-open interval
                           (VL, VU] will be found.
          = 'I': ("Index") the IL-th through IU-th eigenvalues (of the
                           entire matrix) will be found.
[in]N
          N is INTEGER
          The order of the matrix. N > 0.
[in,out]VL
          VL is DOUBLE PRECISION
          If RANGE='V', the lower bound for the eigenvalues.
          Eigenvalues less than or equal to VL, or greater than VU,
          will not be returned.  VL < VU.
          If RANGE='I' or ='A', DLARRE computes bounds on the desired
          part of the spectrum.
[in,out]VU
          VU is DOUBLE PRECISION
          If RANGE='V', the upper bound for the eigenvalues.
          Eigenvalues less than or equal to VL, or greater than VU,
          will not be returned.  VL < VU.
          If RANGE='I' or ='A', DLARRE computes bounds on the desired
          part of the spectrum.
[in]IL
          IL is INTEGER
          If RANGE='I', the index of the
          smallest eigenvalue to be returned.
          1 <= IL <= IU <= N.
[in]IU
          IU is INTEGER
          If RANGE='I', the index of the
          largest eigenvalue to be returned.
          1 <= IL <= IU <= N.
[in,out]D
          D is DOUBLE PRECISION array, dimension (N)
          On entry, the N diagonal elements of the tridiagonal
          matrix T.
          On exit, the N diagonal elements of the diagonal
          matrices D_i.
[in,out]E
          E is DOUBLE PRECISION array, dimension (N)
          On entry, the first (N-1) entries contain the subdiagonal
          elements of the tridiagonal matrix T; E(N) need not be set.
          On exit, E contains the subdiagonal elements of the unit
          bidiagonal matrices L_i. The entries E( ISPLIT( I ) ),
          1 <= I <= NSPLIT, contain the base points sigma_i on output.
[in,out]E2
          E2 is DOUBLE PRECISION array, dimension (N)
          On entry, the first (N-1) entries contain the SQUARES of the
          subdiagonal elements of the tridiagonal matrix T;
          E2(N) need not be set.
          On exit, the entries E2( ISPLIT( I ) ),
          1 <= I <= NSPLIT, have been set to zero
[in]RTOL1
          RTOL1 is DOUBLE PRECISION
[in]RTOL2
          RTOL2 is DOUBLE PRECISION
           Parameters for bisection.
           An interval [LEFT,RIGHT] has converged if
           RIGHT-LEFT < MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
[in]SPLTOL
          SPLTOL is DOUBLE PRECISION
          The threshold for splitting.
[out]NSPLIT
          NSPLIT is INTEGER
          The number of blocks T splits into. 1 <= NSPLIT <= N.
[out]ISPLIT
          ISPLIT is INTEGER array, dimension (N)
          The splitting points, at which T breaks up into blocks.
          The first block consists of rows/columns 1 to ISPLIT(1),
          the second of rows/columns ISPLIT(1)+1 through ISPLIT(2),
          etc., and the NSPLIT-th consists of rows/columns
          ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N.
[out]M
          M is INTEGER
          The total number of eigenvalues (of all L_i D_i L_i^T)
          found.
[out]W
          W is DOUBLE PRECISION array, dimension (N)
          The first M elements contain the eigenvalues. The
          eigenvalues of each of the blocks, L_i D_i L_i^T, are
          sorted in ascending order ( DLARRE may use the
          remaining N-M elements as workspace).
[out]WERR
          WERR is DOUBLE PRECISION array, dimension (N)
          The error bound on the corresponding eigenvalue in W.
[out]WGAP
          WGAP is DOUBLE PRECISION array, dimension (N)
          The separation from the right neighbor eigenvalue in W.
          The gap is only with respect to the eigenvalues of the same block
          as each block has its own representation tree.
          Exception: at the right end of a block we store the left gap
[out]IBLOCK
          IBLOCK is INTEGER array, dimension (N)
          The indices of the blocks (submatrices) associated with the
          corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue
          W(i) belongs to the first block from the top, =2 if W(i)
          belongs to the second block, etc.
[out]INDEXW
          INDEXW is INTEGER array, dimension (N)
          The indices of the eigenvalues within each block (submatrix);
          for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the
          i-th eigenvalue W(i) is the 10-th eigenvalue in block 2
[out]GERS
          GERS is DOUBLE PRECISION array, dimension (2*N)
          The N Gerschgorin intervals (the i-th Gerschgorin interval
          is (GERS(2*i-1), GERS(2*i)).
[out]PIVMIN
          PIVMIN is DOUBLE PRECISION
          The minimum pivot in the Sturm sequence for T.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (6*N)
          Workspace.
[out]IWORK
          IWORK is INTEGER array, dimension (5*N)
          Workspace.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          > 0:  A problem occurred in DLARRE.
          < 0:  One of the called subroutines signaled an internal problem.
                Needs inspection of the corresponding parameter IINFO
                for further information.

          =-1:  Problem in DLARRD.
          = 2:  No base representation could be found in MAXTRY iterations.
                Increasing MAXTRY and recompilation might be a remedy.
          =-3:  Problem in DLARRB when computing the refined root
                representation for DLASQ2.
          =-4:  Problem in DLARRB when preforming bisection on the
                desired part of the spectrum.
          =-5:  Problem in DLASQ2.
          =-6:  Problem in DLASQ2.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  The base representations are required to suffer very little
  element growth and consequently define all their eigenvalues to
  high relative accuracy.
Contributors:
Beresford Parlett, University of California, Berkeley, USA
Jim Demmel, University of California, Berkeley, USA
Inderjit Dhillon, University of Texas, Austin, USA
Osni Marques, LBNL/NERSC, USA
Christof Voemel, University of California, Berkeley, USA

Definition at line 301 of file dlarre.f.

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*
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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 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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
Here is the call graph for this function:
Here is the caller graph for this function: