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

◆ zstemr()

subroutine zstemr ( 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,
complex*16, 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 )

ZSTEMR

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

Purpose:
!>
!> ZSTEMR 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.ZSTEMR 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,
!> ZSTEMR accepts complex workspace to facilitate interoperability
!> with ZUNMTR or ZUPMTR.
!> 
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 COMPLEX*16 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 ZLARRV.
!>                Here, the digit X = ABS( IINFO ) < 10, where IINFO is
!>                the nonzero error code returned by DLARRE or
!>                ZLARRV, 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 334 of file zstemr.f.

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