LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine slarre ( character  RANGE,
integer  N,
real  VL,
real  VU,
integer  IL,
integer  IU,
real, dimension( * )  D,
real, dimension( * )  E,
real, dimension( * )  E2,
real  RTOL1,
real  RTOL2,
real  SPLTOL,
integer  NSPLIT,
integer, dimension( * )  ISPLIT,
integer  M,
real, dimension( * )  W,
real, dimension( * )  WERR,
real, dimension( * )  WGAP,
integer, dimension( * )  IBLOCK,
integer, dimension( * )  INDEXW,
real, dimension( * )  GERS,
real  PIVMIN,
real, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

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

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

Purpose: 
 To find the desired eigenvalues of a given real symmetric
 tridiagonal matrix T, SLARRE 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
 SSTEMR 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 SLASQ2) to
 conpute all and then discard any unwanted one.
 As an added benefit, SLARRE 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 REAL
          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', SLARRE computes bounds on the desired
          part of the spectrum.
[in,out]VU
          VU is REAL
          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', SLARRE 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 REAL 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 REAL 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 REAL 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 REAL
[in]RTOL2
          RTOL2 is REAL
           Parameters for bisection.
           An interval [LEFT,RIGHT] has converged if
           RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
[in]SPLTOL
          SPLTOL is REAL
          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 REAL 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 ( SLARRE may use the
          remaining N-M elements as workspace).
[out]WERR
          WERR is REAL array, dimension (N)
          The error bound on the corresponding eigenvalue in W.
[out]WGAP
          WGAP is REAL 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 REAL 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 REAL
          The minimum pivot in the Sturm sequence for T.
[out]WORK
          WORK is REAL 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 SLARRE.
          < 0:  One of the called subroutines signaled an internal problem.
                Needs inspection of the corresponding parameter IINFO
                for further information.

          =-1:  Problem in SLARRD.
          = 2:  No base representation could be found in MAXTRY iterations.
                Increasing MAXTRY and recompilation might be a remedy.
          =-3:  Problem in SLARRB when computing the refined root
                representation for SLASQ2.
          =-4:  Problem in SLARRB when preforming bisection on the
                desired part of the spectrum.
          =-5:  Problem in SLASQ2.
          =-6:  Problem in SLASQ2.
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 slarre.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: