LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages

◆ sstemr()

subroutine sstemr ( 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,
real, dimension( ldz, * ) z,
integer ldz,
integer nzc,
integer, dimension( * ) isuppz,
logical tryrac,
real, dimension( * ) work,
integer lwork,
integer, dimension( * ) iwork,
integer liwork,
integer info )

SSTEMR

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

Purpose:
!> !> SSTEMR 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: !> Linear Algebra and its Applications, 387(1), pp. 1-28, August 2004. !> - Inderjit Dhillon and Beresford Parlett: SIAM Journal on Matrix Analysis and Applications, Vol. 25, !> 2004. Also LAPACK Working Note 154. !> - Inderjit Dhillon: , !> Computer Science Division Technical Report No. UCB/CSD-97-971, !> UC Berkeley, May 1997. !> !> Further Details !> 1.SSTEMR 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 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 REAL 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 SLARRV. !> Here, the digit X = ABS( IINFO ) < 10, where IINFO is !> the nonzero error code returned by SLARRE or !> SLARRV, 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 317 of file sstemr.f.

320*
321* -- LAPACK computational routine --
322* -- LAPACK is a software package provided by Univ. of Tennessee, --
323* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
324*
325* .. Scalar Arguments ..
326 CHARACTER JOBZ, RANGE
327 LOGICAL TRYRAC
328 INTEGER IL, INFO, IU, LDZ, NZC, LIWORK, LWORK, M, N
329 REAL VL, VU
330* ..
331* .. Array Arguments ..
332 INTEGER ISUPPZ( * ), IWORK( * )
333 REAL D( * ), E( * ), W( * ), WORK( * )
334 REAL Z( LDZ, * )
335* ..
336*
337* =====================================================================
338*
339* .. Parameters ..
340 REAL ZERO, ONE, FOUR, MINRGP
341 parameter( zero = 0.0e0, one = 1.0e0,
342 $ four = 4.0e0,
343 $ minrgp = 3.0e-3 )
344* ..
345* .. Local Scalars ..
346 LOGICAL ALLEIG, INDEIG, LQUERY, VALEIG, WANTZ, ZQUERY,
347 $ LAESWAP
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 REAL 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 REAL SLAMCH, SLANST, SROUNDUP_LWORK
362* ..
363* .. External Subroutines ..
364 EXTERNAL scopy, slae2, slaev2, slarrc, slarre,
365 $ slarrj,
367* ..
368* .. Intrinsic Functions ..
369 INTRINSIC max, min, sqrt
370* ..
371* .. Executable Statements ..
372*
373* Test the input parameters.
374*
375 wantz = lsame( jobz, 'V' )
376 alleig = lsame( range, 'A' )
377 valeig = lsame( range, 'V' )
378 indeig = lsame( range, 'I' )
379*
380 lquery = ( ( lwork.EQ.-1 ).OR.( liwork.EQ.-1 ) )
381 zquery = ( nzc.EQ.-1 )
382 laeswap = .false.
383
384* SSTEMR needs WORK of size 6*N, IWORK of size 3*N.
385* In addition, SLARRE needs WORK of size 6*N, IWORK of size 5*N.
386* Furthermore, SLARRV 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 SLARRE.
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 = slamch( 'Safe minimum' )
438 eps = slamch( '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 ) = sroundup_lwork(lwmin)
446 iwork( 1 ) = liwmin
447*
448 IF( wantz .AND. alleig ) THEN
449 nzcmin = n
450 ELSE IF( wantz .AND. valeig ) THEN
451 CALL slarrc( '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 ) = real( 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( 'SSTEMR', -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 slae2( d(1), e(1), d(2), r1, r2 )
502 ELSE IF( wantz.AND.(.NOT.zquery) ) THEN
503 CALL slaev2( d(1), e(1), d(2), r1, r2, cs, sn )
504 END IF
505* D/S/LAE2 and D/S/LAEV2 outputs satisfy |R1| >= |R2|. However,
506* the following code requires R1 >= R2. Hence, we correct
507* the order of R1, R2, CS, SN if R1 < R2 before further processing.
508 IF( r1.LT.r2 ) THEN
509 e(2) = r1
510 r1 = r2
511 r2 = e(2)
512 laeswap = .true.
513 ENDIF
514 IF( alleig.OR.
515 $ (valeig.AND.(r2.GT.wl).AND.
516 $ (r2.LE.wu)).OR.
517 $ (indeig.AND.(iil.EQ.1)) ) THEN
518 m = m+1
519 w( m ) = r2
520 IF( wantz.AND.(.NOT.zquery) ) THEN
521 IF( laeswap ) THEN
522 z( 1, m ) = cs
523 z( 2, m ) = sn
524 ELSE
525 z( 1, m ) = -sn
526 z( 2, m ) = cs
527 ENDIF
528* Note: At most one of SN and CS can be zero.
529 IF (sn.NE.zero) THEN
530 IF (cs.NE.zero) THEN
531 isuppz(2*m-1) = 1
532 isuppz(2*m) = 2
533 ELSE
534 isuppz(2*m-1) = 1
535 isuppz(2*m) = 1
536 END IF
537 ELSE
538 isuppz(2*m-1) = 2
539 isuppz(2*m) = 2
540 END IF
541 ENDIF
542 ENDIF
543 IF( alleig.OR.
544 $ (valeig.AND.(r1.GT.wl).AND.
545 $ (r1.LE.wu)).OR.
546 $ (indeig.AND.(iiu.EQ.2)) ) THEN
547 m = m+1
548 w( m ) = r1
549 IF( wantz.AND.(.NOT.zquery) ) THEN
550 IF( laeswap ) THEN
551 z( 1, m ) = -sn
552 z( 2, m ) = cs
553 ELSE
554 z( 1, m ) = cs
555 z( 2, m ) = sn
556 ENDIF
557* Note: At most one of SN and CS can be zero.
558 IF (sn.NE.zero) THEN
559 IF (cs.NE.zero) THEN
560 isuppz(2*m-1) = 1
561 isuppz(2*m) = 2
562 ELSE
563 isuppz(2*m-1) = 1
564 isuppz(2*m) = 1
565 END IF
566 ELSE
567 isuppz(2*m-1) = 2
568 isuppz(2*m) = 2
569 END IF
570 ENDIF
571 ENDIF
572 ELSE
573
574* Continue with general N
575
576 indgrs = 1
577 inderr = 2*n + 1
578 indgp = 3*n + 1
579 indd = 4*n + 1
580 inde2 = 5*n + 1
581 indwrk = 6*n + 1
582*
583 iinspl = 1
584 iindbl = n + 1
585 iindw = 2*n + 1
586 iindwk = 3*n + 1
587*
588* Scale matrix to allowable range, if necessary.
589* The allowable range is related to the PIVMIN parameter; see the
590* comments in SLARRD. The preference for scaling small values
591* up is heuristic; we expect users' matrices not to be close to the
592* RMAX threshold.
593*
594 scale = one
595 tnrm = slanst( 'M', n, d, e )
596 IF( tnrm.GT.zero .AND. tnrm.LT.rmin ) THEN
597 scale = rmin / tnrm
598 ELSE IF( tnrm.GT.rmax ) THEN
599 scale = rmax / tnrm
600 END IF
601 IF( scale.NE.one ) THEN
602 CALL sscal( n, scale, d, 1 )
603 CALL sscal( n-1, scale, e, 1 )
604 tnrm = tnrm*scale
605 IF( valeig ) THEN
606* If eigenvalues in interval have to be found,
607* scale (WL, WU] accordingly
608 wl = wl*scale
609 wu = wu*scale
610 ENDIF
611 END IF
612*
613* Compute the desired eigenvalues of the tridiagonal after splitting
614* into smaller subblocks if the corresponding off-diagonal elements
615* are small
616* THRESH is the splitting parameter for SLARRE
617* A negative THRESH forces the old splitting criterion based on the
618* size of the off-diagonal. A positive THRESH switches to splitting
619* which preserves relative accuracy.
620*
621 IF( tryrac ) THEN
622* Test whether the matrix warrants the more expensive relative approach.
623 CALL slarrr( n, d, e, iinfo )
624 ELSE
625* The user does not care about relative accurately eigenvalues
626 iinfo = -1
627 ENDIF
628* Set the splitting criterion
629 IF (iinfo.EQ.0) THEN
630 thresh = eps
631 ELSE
632 thresh = -eps
633* relative accuracy is desired but T does not guarantee it
634 tryrac = .false.
635 ENDIF
636*
637 IF( tryrac ) THEN
638* Copy original diagonal, needed to guarantee relative accuracy
639 CALL scopy(n,d,1,work(indd),1)
640 ENDIF
641* Store the squares of the offdiagonal values of T
642 DO 5 j = 1, n-1
643 work( inde2+j-1 ) = e(j)**2
644 5 CONTINUE
645
646* Set the tolerance parameters for bisection
647 IF( .NOT.wantz ) THEN
648* SLARRE computes the eigenvalues to full precision.
649 rtol1 = four * eps
650 rtol2 = four * eps
651 ELSE
652* SLARRE computes the eigenvalues to less than full precision.
653* SLARRV will refine the eigenvalue approximations, and we can
654* need less accurate initial bisection in SLARRE.
655* Note: these settings do only affect the subset case and SLARRE
656 rtol1 = max( sqrt(eps)*5.0e-2, four * eps )
657 rtol2 = max( sqrt(eps)*5.0e-3, four * eps )
658 ENDIF
659 CALL slarre( range, n, wl, wu, iil, iiu, d, e,
660 $ work(inde2), rtol1, rtol2, thresh, nsplit,
661 $ iwork( iinspl ), m, w, work( inderr ),
662 $ work( indgp ), iwork( iindbl ),
663 $ iwork( iindw ), work( indgrs ), pivmin,
664 $ work( indwrk ), iwork( iindwk ), iinfo )
665 IF( iinfo.NE.0 ) THEN
666 info = 10 + abs( iinfo )
667 RETURN
668 END IF
669* Note that if RANGE .NE. 'V', SLARRE computes bounds on the desired
670* part of the spectrum. All desired eigenvalues are contained in
671* (WL,WU]
672
673
674 IF( wantz ) THEN
675*
676* Compute the desired eigenvectors corresponding to the computed
677* eigenvalues
678*
679 CALL slarrv( n, wl, wu, d, e,
680 $ pivmin, iwork( iinspl ), m,
681 $ 1, m, minrgp, rtol1, rtol2,
682 $ w, work( inderr ), work( indgp ), iwork( iindbl ),
683 $ iwork( iindw ), work( indgrs ), z, ldz,
684 $ isuppz, work( indwrk ), iwork( iindwk ), iinfo )
685 IF( iinfo.NE.0 ) THEN
686 info = 20 + abs( iinfo )
687 RETURN
688 END IF
689 ELSE
690* SLARRE computes eigenvalues of the (shifted) root representation
691* SLARRV returns the eigenvalues of the unshifted matrix.
692* However, if the eigenvectors are not desired by the user, we need
693* to apply the corresponding shifts from SLARRE to obtain the
694* eigenvalues of the original matrix.
695 DO 20 j = 1, m
696 itmp = iwork( iindbl+j-1 )
697 w( j ) = w( j ) + e( iwork( iinspl+itmp-1 ) )
698 20 CONTINUE
699 END IF
700*
701
702 IF ( tryrac ) THEN
703* Refine computed eigenvalues so that they are relatively accurate
704* with respect to the original matrix T.
705 ibegin = 1
706 wbegin = 1
707 DO 39 jblk = 1, iwork( iindbl+m-1 )
708 iend = iwork( iinspl+jblk-1 )
709 in = iend - ibegin + 1
710 wend = wbegin - 1
711* check if any eigenvalues have to be refined in this block
712 36 CONTINUE
713 IF( wend.LT.m ) THEN
714 IF( iwork( iindbl+wend ).EQ.jblk ) THEN
715 wend = wend + 1
716 GO TO 36
717 END IF
718 END IF
719 IF( wend.LT.wbegin ) THEN
720 ibegin = iend + 1
721 GO TO 39
722 END IF
723
724 offset = iwork(iindw+wbegin-1)-1
725 ifirst = iwork(iindw+wbegin-1)
726 ilast = iwork(iindw+wend-1)
727 rtol2 = four * eps
728 CALL slarrj( in,
729 $ work(indd+ibegin-1), work(inde2+ibegin-1),
730 $ ifirst, ilast, rtol2, offset, w(wbegin),
731 $ work( inderr+wbegin-1 ),
732 $ work( indwrk ), iwork( iindwk ), pivmin,
733 $ tnrm, iinfo )
734 ibegin = iend + 1
735 wbegin = wend + 1
736 39 CONTINUE
737 ENDIF
738*
739* If matrix was scaled, then rescale eigenvalues appropriately.
740*
741 IF( scale.NE.one ) THEN
742 CALL sscal( m, one / scale, w, 1 )
743 END IF
744 END IF
745*
746* If eigenvalues are not in increasing order, then sort them,
747* possibly along with eigenvectors.
748*
749 IF( nsplit.GT.1 .OR. n.EQ.2 ) THEN
750 IF( .NOT. wantz ) THEN
751 CALL slasrt( 'I', m, w, iinfo )
752 IF( iinfo.NE.0 ) THEN
753 info = 3
754 RETURN
755 END IF
756 ELSE
757 DO 60 j = 1, m - 1
758 i = 0
759 tmp = w( j )
760 DO 50 jj = j + 1, m
761 IF( w( jj ).LT.tmp ) THEN
762 i = jj
763 tmp = w( jj )
764 END IF
765 50 CONTINUE
766 IF( i.NE.0 ) THEN
767 w( i ) = w( j )
768 w( j ) = tmp
769 IF( wantz ) THEN
770 CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
771 itmp = isuppz( 2*i-1 )
772 isuppz( 2*i-1 ) = isuppz( 2*j-1 )
773 isuppz( 2*j-1 ) = itmp
774 itmp = isuppz( 2*i )
775 isuppz( 2*i ) = isuppz( 2*j )
776 isuppz( 2*j ) = itmp
777 END IF
778 END IF
779 60 CONTINUE
780 END IF
781 ENDIF
782*
783*
784 work( 1 ) = sroundup_lwork(lwmin)
785 iwork( 1 ) = liwmin
786 RETURN
787*
788* End of SSTEMR
789*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine slae2(a, b, c, rt1, rt2)
SLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix.
Definition slae2.f:100
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:118
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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:98
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:135
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:303
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:166
subroutine slarrr(n, d, e, info)
SLARRR performs tests to decide whether the symmetric tridiagonal matrix T warrants expensive computa...
Definition slarrr.f:92
subroutine slarrv(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)
SLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenvalues ...
Definition slarrv.f:290
subroutine slasrt(id, n, d, info)
SLASRT sorts numbers in increasing or decreasing order.
Definition slasrt.f:86
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: