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

◆ dstemr()

subroutine dstemr ( character  jobz,
character  range,
integer  n,
double precision, dimension( * )  d,
double precision, dimension( * )  e,
double precision  vl,
double precision  vu,
integer  il,
integer  iu,
integer  m,
double precision, dimension( * )  w,
double precision, dimension( ldz, * )  z,
integer  ldz,
integer  nzc,
integer, dimension( * )  isuppz,
logical  tryrac,
double precision, dimension( * )  work,
integer  lwork,
integer, dimension( * )  iwork,
integer  liwork,
integer  info 
)

DSTEMR

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

Purpose:
 DSTEMR computes selected eigenvalues and, optionally, eigenvectors
 of a real symmetric tridiagonal matrix T. Any such unreduced matrix has
 a well defined set of pairwise different real eigenvalues, the corresponding
 real eigenvectors are pairwise orthogonal.

 The spectrum may be computed either completely or partially by specifying
 either an interval (VL,VU] or a range of indices IL:IU for the desired
 eigenvalues.

 Depending on the number of desired eigenvalues, these are computed either
 by bisection or the dqds algorithm. Numerically orthogonal eigenvectors are
 computed by the use of various suitable L D L^T factorizations near clusters
 of close eigenvalues (referred to as RRRs, Relatively Robust
 Representations). An informal sketch of the algorithm follows.

 For each unreduced block (submatrix) of T,
    (a) Compute T - sigma I  = L D L^T, so that L and D
        define all the wanted eigenvalues to high relative accuracy.
        This means that small relative changes in the entries of D and L
        cause only small relative changes in the eigenvalues and
        eigenvectors. The standard (unfactored) representation of the
        tridiagonal matrix T does not have this property in general.
    (b) Compute the eigenvalues to suitable accuracy.
        If the eigenvectors are desired, the algorithm attains full
        accuracy of the computed eigenvalues only right before
        the corresponding vectors have to be computed, see steps c) and d).
    (c) For each cluster of close eigenvalues, select a new
        shift close to the cluster, find a new factorization, and refine
        the shifted eigenvalues to suitable accuracy.
    (d) For each eigenvalue with a large enough relative separation compute
        the corresponding eigenvector by forming a rank revealing twisted
        factorization. Go back to (c) for any clusters that remain.

 For more details, see:
 - Inderjit S. Dhillon and Beresford N. Parlett: "Multiple representations
   to compute orthogonal eigenvectors of symmetric tridiagonal matrices,"
   Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004.
 - Inderjit Dhillon and Beresford Parlett: "Orthogonal Eigenvectors and
   Relative Gaps," SIAM Journal on Matrix Analysis and Applications, Vol. 25,
   2004.  Also LAPACK Working Note 154.
 - Inderjit Dhillon: "A new O(n^2) algorithm for the symmetric
   tridiagonal eigenvalue/eigenvector problem",
   Computer Science Division Technical Report No. UCB/CSD-97-971,
   UC Berkeley, May 1997.

 Further Details
 1.DSTEMR works only on machines which follow IEEE-754
 floating-point standard in their handling of infinities and NaNs.
 This permits the use of efficient inner loops avoiding a check for
 zero divisors.
Parameters
[in]JOBZ
          JOBZ is CHARACTER*1
          = 'N':  Compute eigenvalues only;
          = 'V':  Compute eigenvalues and eigenvectors.
[in]RANGE
          RANGE is CHARACTER*1
          = 'A': all eigenvalues will be found.
          = 'V': all eigenvalues in the half-open interval (VL,VU]
                 will be found.
          = 'I': the IL-th through IU-th eigenvalues will be found.
[in]N
          N is INTEGER
          The order of the matrix.  N >= 0.
[in,out]D
          D is DOUBLE PRECISION array, dimension (N)
          On entry, the N diagonal elements of the tridiagonal matrix
          T. On exit, D is overwritten.
[in,out]E
          E is DOUBLE PRECISION array, dimension (N)
          On entry, the (N-1) subdiagonal elements of the tridiagonal
          matrix T in elements 1 to N-1 of E. E(N) need not be set on
          input, but is used internally as workspace.
          On exit, E is overwritten.
[in]VL
          VL is DOUBLE PRECISION

          If RANGE='V', the lower bound of the interval to
          be searched for eigenvalues. VL < VU.
          Not referenced if RANGE = 'A' or 'I'.
[in]VU
          VU is DOUBLE PRECISION

          If RANGE='V', the upper bound of the interval to
          be searched for eigenvalues. VL < VU.
          Not referenced if RANGE = 'A' or 'I'.
[in]IL
          IL is INTEGER

          If RANGE='I', the index of the
          smallest eigenvalue to be returned.
          1 <= IL <= IU <= N, if N > 0.
          Not referenced if RANGE = 'A' or 'V'.
[in]IU
          IU is INTEGER

          If RANGE='I', the index of the
          largest eigenvalue to be returned.
          1 <= IL <= IU <= N, if N > 0.
          Not referenced if RANGE = 'A' or 'V'.
[out]M
          M is INTEGER
          The total number of eigenvalues found.  0 <= M <= N.
          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
[out]W
          W is DOUBLE PRECISION array, dimension (N)
          The first M elements contain the selected eigenvalues in
          ascending order.
[out]Z
          Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M) )
          If JOBZ = 'V', and if INFO = 0, then the first M columns of Z
          contain the orthonormal eigenvectors of the matrix T
          corresponding to the selected eigenvalues, with the i-th
          column of Z holding the eigenvector associated with W(i).
          If JOBZ = 'N', then Z is not referenced.
          Note: the user must ensure that at least max(1,M) columns are
          supplied in the array Z; if RANGE = 'V', the exact value of M
          is not known in advance and can be computed with a workspace
          query by setting NZC = -1, see below.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z.  LDZ >= 1, and if
          JOBZ = 'V', then LDZ >= max(1,N).
[in]NZC
          NZC is INTEGER
          The number of eigenvectors to be held in the array Z.
          If RANGE = 'A', then NZC >= max(1,N).
          If RANGE = 'V', then NZC >= the number of eigenvalues in (VL,VU].
          If RANGE = 'I', then NZC >= IU-IL+1.
          If NZC = -1, then a workspace query is assumed; the
          routine calculates the number of columns of the array Z that
          are needed to hold the eigenvectors.
          This value is returned as the first entry of the Z array, and
          no error message related to NZC is issued by XERBLA.
[out]ISUPPZ
          ISUPPZ is INTEGER array, dimension ( 2*max(1,M) )
          The support of the eigenvectors in Z, i.e., the indices
          indicating the nonzero elements in Z. The i-th computed eigenvector
          is nonzero only in elements ISUPPZ( 2*i-1 ) through
          ISUPPZ( 2*i ). This is relevant in the case when the matrix
          is split. ISUPPZ is only accessed when JOBZ is 'V' and N > 0.
[in,out]TRYRAC
          TRYRAC is LOGICAL
          If TRYRAC = .TRUE., indicates that the code should check whether
          the tridiagonal matrix defines its eigenvalues to high relative
          accuracy.  If so, the code uses relative-accuracy preserving
          algorithms that might be (a bit) slower depending on the matrix.
          If the matrix does not define its eigenvalues to high relative
          accuracy, the code can uses possibly faster algorithms.
          If TRYRAC = .FALSE., the code is not required to guarantee
          relatively accurate eigenvalues and can use the fastest possible
          techniques.
          On exit, a .TRUE. TRYRAC will be set to .FALSE. if the matrix
          does not define its eigenvalues to high relative accuracy.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (LWORK)
          On exit, if INFO = 0, WORK(1) returns the optimal
          (and minimal) LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= max(1,18*N)
          if JOBZ = 'V', and LWORK >= max(1,12*N) if JOBZ = 'N'.
          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]IWORK
          IWORK is INTEGER array, dimension (LIWORK)
          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
[in]LIWORK
          LIWORK is INTEGER
          The dimension of the array IWORK.  LIWORK >= max(1,10*N)
          if the eigenvectors are desired, and LIWORK >= max(1,8*N)
          if only the eigenvalues are to be computed.
          If LIWORK = -1, then a workspace query is assumed; the
          routine only calculates the optimal size of the IWORK array,
          returns this value as the first entry of the IWORK array, and
          no error message related to LIWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          On exit, INFO
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = 1X, internal error in DLARRE,
                if INFO = 2X, internal error in DLARRV.
                Here, the digit X = ABS( IINFO ) < 10, where IINFO is
                the nonzero error code returned by DLARRE or
                DLARRV, respectively.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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
Aravindh Krishnamoorthy, FAU, Erlangen, Germany

Definition at line 319 of file dstemr.f.

322*
323* -- LAPACK computational routine --
324* -- LAPACK is a software package provided by Univ. of Tennessee, --
325* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
326*
327* .. Scalar Arguments ..
328 CHARACTER JOBZ, RANGE
329 LOGICAL TRYRAC
330 INTEGER IL, INFO, IU, LDZ, NZC, LIWORK, LWORK, M, N
331 DOUBLE PRECISION VL, VU
332* ..
333* .. Array Arguments ..
334 INTEGER ISUPPZ( * ), IWORK( * )
335 DOUBLE PRECISION D( * ), E( * ), W( * ), WORK( * )
336 DOUBLE PRECISION Z( LDZ, * )
337* ..
338*
339* =====================================================================
340*
341* .. Parameters ..
342 DOUBLE PRECISION ZERO, ONE, FOUR, MINRGP
343 parameter( zero = 0.0d0, one = 1.0d0,
344 $ four = 4.0d0,
345 $ minrgp = 1.0d-3 )
346* ..
347* .. Local Scalars ..
348 LOGICAL ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY,
349 $ LAESWAP
350 INTEGER I, IBEGIN, IEND, IFIRST, IIL, IINDBL, IINDW,
351 $ IINDWK, IINFO, IINSPL, IIU, ILAST, IN, INDD,
352 $ INDE2, INDERR, INDGP, INDGRS, INDWRK, ITMP,
353 $ ITMP2, J, JBLK, JJ, LIWMIN, LWMIN, NSPLIT,
354 $ NZCMIN, OFFSET, WBEGIN, WEND
355 DOUBLE PRECISION BIGNUM, CS, EPS, PIVMIN, R1, R2, RMAX, RMIN,
356 $ RTOL1, RTOL2, SAFMIN, SCALE, SMLNUM, SN,
357 $ THRESH, TMP, TNRM, WL, WU
358* ..
359* ..
360* .. External Functions ..
361 LOGICAL LSAME
362 DOUBLE PRECISION DLAMCH, DLANST
363 EXTERNAL lsame, dlamch, dlanst
364* ..
365* .. External Subroutines ..
366 EXTERNAL dcopy, dlae2, dlaev2, dlarrc, dlarre, dlarrj,
368* ..
369* .. Intrinsic Functions ..
370 INTRINSIC max, min, sqrt
371
372
373* ..
374* .. Executable Statements ..
375*
376* Test the input parameters.
377*
378 wantz = lsame( jobz, 'V' )
379 alleig = lsame( range, 'A' )
380 valeig = lsame( range, 'V' )
381 indeig = lsame( range, 'I' )
382*
383 lquery = ( ( lwork.EQ.-1 ).OR.( liwork.EQ.-1 ) )
384 zquery = ( nzc.EQ.-1 )
385 laeswap = .false.
386
387* DSTEMR needs WORK of size 6*N, IWORK of size 3*N.
388* In addition, DLARRE needs WORK of size 6*N, IWORK of size 5*N.
389* Furthermore, DLARRV needs WORK of size 12*N, IWORK of size 7*N.
390 IF( wantz ) THEN
391 lwmin = 18*n
392 liwmin = 10*n
393 ELSE
394* need less workspace if only the eigenvalues are wanted
395 lwmin = 12*n
396 liwmin = 8*n
397 ENDIF
398
399 wl = zero
400 wu = zero
401 iil = 0
402 iiu = 0
403 nsplit = 0
404
405 IF( valeig ) THEN
406* We do not reference VL, VU in the cases RANGE = 'I','A'
407* The interval (WL, WU] contains all the wanted eigenvalues.
408* It is either given by the user or computed in DLARRE.
409 wl = vl
410 wu = vu
411 ELSEIF( indeig ) THEN
412* We do not reference IL, IU in the cases RANGE = 'V','A'
413 iil = il
414 iiu = iu
415 ENDIF
416*
417 info = 0
418 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
419 info = -1
420 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
421 info = -2
422 ELSE IF( n.LT.0 ) THEN
423 info = -3
424 ELSE IF( valeig .AND. n.GT.0 .AND. wu.LE.wl ) THEN
425 info = -7
426 ELSE IF( indeig .AND. ( iil.LT.1 .OR. iil.GT.n ) ) THEN
427 info = -8
428 ELSE IF( indeig .AND. ( iiu.LT.iil .OR. iiu.GT.n ) ) THEN
429 info = -9
430 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
431 info = -13
432 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
433 info = -17
434 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
435 info = -19
436 END IF
437*
438* Get machine constants.
439*
440 safmin = dlamch( 'Safe minimum' )
441 eps = dlamch( 'Precision' )
442 smlnum = safmin / eps
443 bignum = one / smlnum
444 rmin = sqrt( smlnum )
445 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
446*
447 IF( info.EQ.0 ) THEN
448 work( 1 ) = lwmin
449 iwork( 1 ) = liwmin
450*
451 IF( wantz .AND. alleig ) THEN
452 nzcmin = n
453 ELSE IF( wantz .AND. valeig ) THEN
454 CALL dlarrc( 'T', n, vl, vu, d, e, safmin,
455 $ nzcmin, itmp, itmp2, info )
456 ELSE IF( wantz .AND. indeig ) THEN
457 nzcmin = iiu-iil+1
458 ELSE
459* WANTZ .EQ. FALSE.
460 nzcmin = 0
461 ENDIF
462 IF( zquery .AND. info.EQ.0 ) THEN
463 z( 1,1 ) = nzcmin
464 ELSE IF( nzc.LT.nzcmin .AND. .NOT.zquery ) THEN
465 info = -14
466 END IF
467 END IF
468
469 IF( info.NE.0 ) THEN
470*
471 CALL xerbla( 'DSTEMR', -info )
472*
473 RETURN
474 ELSE IF( lquery .OR. zquery ) THEN
475 RETURN
476 END IF
477*
478* Handle N = 0, 1, and 2 cases immediately
479*
480 m = 0
481 IF( n.EQ.0 )
482 $ RETURN
483*
484 IF( n.EQ.1 ) THEN
485 IF( alleig .OR. indeig ) THEN
486 m = 1
487 w( 1 ) = d( 1 )
488 ELSE
489 IF( wl.LT.d( 1 ) .AND. wu.GE.d( 1 ) ) THEN
490 m = 1
491 w( 1 ) = d( 1 )
492 END IF
493 END IF
494 IF( wantz.AND.(.NOT.zquery) ) THEN
495 z( 1, 1 ) = one
496 isuppz(1) = 1
497 isuppz(2) = 1
498 END IF
499 RETURN
500 END IF
501*
502 IF( n.EQ.2 ) THEN
503 IF( .NOT.wantz ) THEN
504 CALL dlae2( d(1), e(1), d(2), r1, r2 )
505 ELSE IF( wantz.AND.(.NOT.zquery) ) THEN
506 CALL dlaev2( d(1), e(1), d(2), r1, r2, cs, sn )
507 END IF
508* D/S/LAE2 and D/S/LAEV2 outputs satisfy |R1| >= |R2|. However,
509* the following code requires R1 >= R2. Hence, we correct
510* the order of R1, R2, CS, SN if R1 < R2 before further processing.
511 IF( r1.LT.r2 ) THEN
512 e(2) = r1
513 r1 = r2
514 r2 = e(2)
515 laeswap = .true.
516 ENDIF
517 IF( alleig.OR.
518 $ (valeig.AND.(r2.GT.wl).AND.
519 $ (r2.LE.wu)).OR.
520 $ (indeig.AND.(iil.EQ.1)) ) THEN
521 m = m+1
522 w( m ) = r2
523 IF( wantz.AND.(.NOT.zquery) ) THEN
524 IF( laeswap ) THEN
525 z( 1, m ) = cs
526 z( 2, m ) = sn
527 ELSE
528 z( 1, m ) = -sn
529 z( 2, m ) = cs
530 ENDIF
531* Note: At most one of SN and CS can be zero.
532 IF (sn.NE.zero) THEN
533 IF (cs.NE.zero) THEN
534 isuppz(2*m-1) = 1
535 isuppz(2*m) = 2
536 ELSE
537 isuppz(2*m-1) = 1
538 isuppz(2*m) = 1
539 END IF
540 ELSE
541 isuppz(2*m-1) = 2
542 isuppz(2*m) = 2
543 END IF
544 ENDIF
545 ENDIF
546 IF( alleig.OR.
547 $ (valeig.AND.(r1.GT.wl).AND.
548 $ (r1.LE.wu)).OR.
549 $ (indeig.AND.(iiu.EQ.2)) ) THEN
550 m = m+1
551 w( m ) = r1
552 IF( wantz.AND.(.NOT.zquery) ) THEN
553 IF( laeswap ) THEN
554 z( 1, m ) = -sn
555 z( 2, m ) = cs
556 ELSE
557 z( 1, m ) = cs
558 z( 2, m ) = sn
559 ENDIF
560* Note: At most one of SN and CS can be zero.
561 IF (sn.NE.zero) THEN
562 IF (cs.NE.zero) THEN
563 isuppz(2*m-1) = 1
564 isuppz(2*m) = 2
565 ELSE
566 isuppz(2*m-1) = 1
567 isuppz(2*m) = 1
568 END IF
569 ELSE
570 isuppz(2*m-1) = 2
571 isuppz(2*m) = 2
572 END IF
573 ENDIF
574 ENDIF
575
576 ELSE
577
578* Continue with general N
579
580 indgrs = 1
581 inderr = 2*n + 1
582 indgp = 3*n + 1
583 indd = 4*n + 1
584 inde2 = 5*n + 1
585 indwrk = 6*n + 1
586*
587 iinspl = 1
588 iindbl = n + 1
589 iindw = 2*n + 1
590 iindwk = 3*n + 1
591*
592* Scale matrix to allowable range, if necessary.
593* The allowable range is related to the PIVMIN parameter; see the
594* comments in DLARRD. The preference for scaling small values
595* up is heuristic; we expect users' matrices not to be close to the
596* RMAX threshold.
597*
598 scale = one
599 tnrm = dlanst( 'M', n, d, e )
600 IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
601 scale = rmin / tnrm
602 ELSE IF( tnrm.GT.rmax ) THEN
603 scale = rmax / tnrm
604 END IF
605 IF( scale.NE.one ) THEN
606 CALL dscal( n, scale, d, 1 )
607 CALL dscal( n-1, scale, e, 1 )
608 tnrm = tnrm*scale
609 IF( valeig ) THEN
610* If eigenvalues in interval have to be found,
611* scale (WL, WU] accordingly
612 wl = wl*scale
613 wu = wu*scale
614 ENDIF
615 END IF
616*
617* Compute the desired eigenvalues of the tridiagonal after splitting
618* into smaller subblocks if the corresponding off-diagonal elements
619* are small
620* THRESH is the splitting parameter for DLARRE
621* A negative THRESH forces the old splitting criterion based on the
622* size of the off-diagonal. A positive THRESH switches to splitting
623* which preserves relative accuracy.
624*
625 IF( tryrac ) THEN
626* Test whether the matrix warrants the more expensive relative approach.
627 CALL dlarrr( n, d, e, iinfo )
628 ELSE
629* The user does not care about relative accurately eigenvalues
630 iinfo = -1
631 ENDIF
632* Set the splitting criterion
633 IF (iinfo.EQ.0) THEN
634 thresh = eps
635 ELSE
636 thresh = -eps
637* relative accuracy is desired but T does not guarantee it
638 tryrac = .false.
639 ENDIF
640*
641 IF( tryrac ) THEN
642* Copy original diagonal, needed to guarantee relative accuracy
643 CALL dcopy(n,d,1,work(indd),1)
644 ENDIF
645* Store the squares of the offdiagonal values of T
646 DO 5 j = 1, n-1
647 work( inde2+j-1 ) = e(j)**2
648 5 CONTINUE
649
650* Set the tolerance parameters for bisection
651 IF( .NOT.wantz ) THEN
652* DLARRE computes the eigenvalues to full precision.
653 rtol1 = four * eps
654 rtol2 = four * eps
655 ELSE
656* DLARRE computes the eigenvalues to less than full precision.
657* DLARRV will refine the eigenvalue approximations, and we can
658* need less accurate initial bisection in DLARRE.
659* Note: these settings do only affect the subset case and DLARRE
660 rtol1 = sqrt(eps)
661 rtol2 = max( sqrt(eps)*5.0d-3, four * eps )
662 ENDIF
663 CALL dlarre( range, n, wl, wu, iil, iiu, d, e,
664 $ work(inde2), rtol1, rtol2, thresh, nsplit,
665 $ iwork( iinspl ), m, w, work( inderr ),
666 $ work( indgp ), iwork( iindbl ),
667 $ iwork( iindw ), work( indgrs ), pivmin,
668 $ work( indwrk ), iwork( iindwk ), iinfo )
669 IF( iinfo.NE.0 ) THEN
670 info = 10 + abs( iinfo )
671 RETURN
672 END IF
673* Note that if RANGE .NE. 'V', DLARRE computes bounds on the desired
674* part of the spectrum. All desired eigenvalues are contained in
675* (WL,WU]
676
677
678 IF( wantz ) THEN
679*
680* Compute the desired eigenvectors corresponding to the computed
681* eigenvalues
682*
683 CALL dlarrv( n, wl, wu, d, e,
684 $ pivmin, iwork( iinspl ), m,
685 $ 1, m, minrgp, rtol1, rtol2,
686 $ w, work( inderr ), work( indgp ), iwork( iindbl ),
687 $ iwork( iindw ), work( indgrs ), z, ldz,
688 $ isuppz, work( indwrk ), iwork( iindwk ), iinfo )
689 IF( iinfo.NE.0 ) THEN
690 info = 20 + abs( iinfo )
691 RETURN
692 END IF
693 ELSE
694* DLARRE computes eigenvalues of the (shifted) root representation
695* DLARRV returns the eigenvalues of the unshifted matrix.
696* However, if the eigenvectors are not desired by the user, we need
697* to apply the corresponding shifts from DLARRE to obtain the
698* eigenvalues of the original matrix.
699 DO 20 j = 1, m
700 itmp = iwork( iindbl+j-1 )
701 w( j ) = w( j ) + e( iwork( iinspl+itmp-1 ) )
702 20 CONTINUE
703 END IF
704*
705
706 IF ( tryrac ) THEN
707* Refine computed eigenvalues so that they are relatively accurate
708* with respect to the original matrix T.
709 ibegin = 1
710 wbegin = 1
711 DO 39 jblk = 1, iwork( iindbl+m-1 )
712 iend = iwork( iinspl+jblk-1 )
713 in = iend - ibegin + 1
714 wend = wbegin - 1
715* check if any eigenvalues have to be refined in this block
716 36 CONTINUE
717 IF( wend.LT.m ) THEN
718 IF( iwork( iindbl+wend ).EQ.jblk ) THEN
719 wend = wend + 1
720 GO TO 36
721 END IF
722 END IF
723 IF( wend.LT.wbegin ) THEN
724 ibegin = iend + 1
725 GO TO 39
726 END IF
727
728 offset = iwork(iindw+wbegin-1)-1
729 ifirst = iwork(iindw+wbegin-1)
730 ilast = iwork(iindw+wend-1)
731 rtol2 = four * eps
732 CALL dlarrj( in,
733 $ work(indd+ibegin-1), work(inde2+ibegin-1),
734 $ ifirst, ilast, rtol2, offset, w(wbegin),
735 $ work( inderr+wbegin-1 ),
736 $ work( indwrk ), iwork( iindwk ), pivmin,
737 $ tnrm, iinfo )
738 ibegin = iend + 1
739 wbegin = wend + 1
740 39 CONTINUE
741 ENDIF
742*
743* If matrix was scaled, then rescale eigenvalues appropriately.
744*
745 IF( scale.NE.one ) THEN
746 CALL dscal( m, one / scale, w, 1 )
747 END IF
748
749 END IF
750
751*
752* If eigenvalues are not in increasing order, then sort them,
753* possibly along with eigenvectors.
754*
755 IF( nsplit.GT.1 .OR. n.EQ.2 ) THEN
756 IF( .NOT. wantz ) THEN
757 CALL dlasrt( 'I', m, w, iinfo )
758 IF( iinfo.NE.0 ) THEN
759 info = 3
760 RETURN
761 END IF
762 ELSE
763 DO 60 j = 1, m - 1
764 i = 0
765 tmp = w( j )
766 DO 50 jj = j + 1, m
767 IF( w( jj ).LT.tmp ) THEN
768 i = jj
769 tmp = w( jj )
770 END IF
771 50 CONTINUE
772 IF( i.NE.0 ) THEN
773 w( i ) = w( j )
774 w( j ) = tmp
775 IF( wantz ) THEN
776 CALL dswap( n, z( 1, i ), 1, z( 1, j ), 1 )
777 itmp = isuppz( 2*i-1 )
778 isuppz( 2*i-1 ) = isuppz( 2*j-1 )
779 isuppz( 2*j-1 ) = itmp
780 itmp = isuppz( 2*i )
781 isuppz( 2*i ) = isuppz( 2*j )
782 isuppz( 2*j ) = itmp
783 END IF
784 END IF
785 60 CONTINUE
786 END IF
787 ENDIF
788*
789*
790 work( 1 ) = lwmin
791 iwork( 1 ) = liwmin
792 RETURN
793*
794* End of DSTEMR
795*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dlae2(a, b, c, rt1, rt2)
DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition dlae2.f:102
subroutine dlaev2(a, b, c, rt1, rt2, cs1, sn1)
DLAEV2 computes the eigenvalues and eigenvectors of a 2-by-2 symmetric/Hermitian matrix.
Definition dlaev2.f:120
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
double precision function dlanst(norm, n, d, e)
DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition dlanst.f:100
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 dlarre(range, n, vl, vu, il, iu, d, e, e2, rtol1, rtol2, spltol, nsplit, isplit, m, w, werr, wgap, iblock, indexw, gers, pivmin, work, iwork, info)
DLARRE given the tridiagonal matrix T, sets small off-diagonal elements to zero and for each unreduce...
Definition dlarre.f:305
subroutine dlarrj(n, d, e2, ifirst, ilast, rtol, offset, w, werr, work, iwork, pivmin, spdiam, info)
DLARRJ performs refinement of the initial estimates of the eigenvalues of the matrix T.
Definition dlarrj.f:168
subroutine dlarrr(n, d, e, info)
DLARRR performs tests to decide whether the symmetric tridiagonal matrix T warrants expensive computa...
Definition dlarrr.f:94
subroutine dlarrv(n, vl, vu, d, l, pivmin, isplit, m, dol, dou, minrgp, rtol1, rtol2, w, werr, wgap, iblock, indexw, gers, z, ldz, isuppz, work, iwork, info)
DLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues ...
Definition dlarrv.f:292
subroutine dlasrt(id, n, d, info)
DLASRT sorts numbers in increasing or decreasing order.
Definition dlasrt.f:88
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: