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

◆ cstemr()

subroutine cstemr ( character  JOBZ,
character  RANGE,
integer  N,
real, dimension( * )  D,
real, dimension( * )  E,
real  VL,
real  VU,
integer  IL,
integer  IU,
integer  M,
real, dimension( * )  W,
complex, dimension( ldz, * )  Z,
integer  LDZ,
integer  NZC,
integer, dimension( * )  ISUPPZ,
logical  TRYRAC,
real, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
integer  LIWORK,
integer  INFO 
)

CSTEMR

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

Purpose:
 CSTEMR 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.CSTEMR 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.

 2. LAPACK routines can be used to reduce a complex Hermitean matrix to
 real symmetric tridiagonal form.

 (Any complex Hermitean tridiagonal matrix has real values on its diagonal
 and potentially complex numbers on its off-diagonals. By applying a
 similarity transform with an appropriate diagonal matrix
 diag(1,e^{i \phy_1}, ... , e^{i \phy_{n-1}}), the complex Hermitean
 matrix can be transformed into a real symmetric matrix and complex
 arithmetic can be entirely avoided.)

 While the eigenvectors of the real symmetric tridiagonal matrix are real,
 the eigenvectors of original complex Hermitean matrix have complex entries
 in general.
 Since LAPACK drivers overwrite the matrix data with the eigenvectors,
 CSTEMR accepts complex workspace to facilitate interoperability
 with CUNMTR or CUPMTR.
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 REAL array, dimension (N)
          On entry, the N diagonal elements of the tridiagonal matrix
          T. On exit, D is overwritten.
[in,out]E
          E is REAL 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 REAL

          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 REAL

          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 REAL array, dimension (N)
          The first M elements contain the selected eigenvalues in
          ascending order.
[out]Z
          Z is COMPLEX 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 REAL 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 SLARRE,
                if INFO = 2X, internal error in CLARRV.
                Here, the digit X = ABS( IINFO ) < 10, where IINFO is
                the nonzero error code returned by SLARRE or
                CLARRV, 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 335 of file cstemr.f.

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