LAPACK 3.11.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

Definition at line 318 of file dstemr.f.

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