LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ 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
 conpute 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.
Date
June 2016
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 307 of file dlarre.f.

307 *
308 * -- LAPACK auxiliary routine (version 3.8.0) --
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, dlarrk
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 *
Here is the call graph for this function:
Here is the caller graph for this function:
dlarrd
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
dlarrk
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
dlarrb
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
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:84
dlasq2
subroutine dlasq2(N, Z, INFO)
DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated ...
Definition: dlasq2.f:114
dlarra
subroutine dlarra(N, D, E, E2, SPLTOL, TNRM, NSPLIT, ISPLIT, INFO)
DLARRA computes the splitting points with the specified threshold.
Definition: dlarra.f:138
dlarnv
subroutine dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: dlarnv.f:99
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
dlarrc
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
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:70