LAPACK  3.10.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.
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 *
371 * Quick return if possible
372 *
373  IF( n.LE.0 ) THEN
374  RETURN
375  END IF
376 *
377 * Decode RANGE
378 *
379  IF( lsame( range, 'A' ) ) THEN
380  irange = allrng
381  ELSE IF( lsame( range, 'V' ) ) THEN
382  irange = valrng
383  ELSE IF( lsame( range, 'I' ) ) THEN
384  irange = indrng
385  END IF
386 
387  m = 0
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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
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 dlarnv(IDIST, ISEED, N, X)
DLARNV returns a vector of random numbers from a uniform or normal distribution.
Definition: dlarnv.f:97
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 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 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
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dlasq2(N, Z, INFO)
DLASQ2 computes all the eigenvalues of the symmetric positive definite tridiagonal matrix associated ...
Definition: dlasq2.f:112
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
Here is the call graph for this function:
Here is the caller graph for this function: