311 SUBROUTINE dstemr( JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
312 $ m, w, z, ldz, nzc, isuppz, tryrac, work, lwork,
313 $ iwork, liwork, info )
321 CHARACTER jobz, range
323 INTEGER il, info, iu, ldz, nzc, liwork, lwork, m, n
324 DOUBLE PRECISION vl, vu
327 INTEGER isuppz( * ), iwork( * )
328 DOUBLE PRECISION d( * ), e( * ), w( * ), work( * )
329 DOUBLE PRECISION z( ldz, * )
335 DOUBLE PRECISION zero, one, four, minrgp
336 parameter( zero = 0.0d0, one = 1.0d0,
341 LOGICAL alleig, indeig, lquery, valeig, wantz, zquery
342 INTEGER i, ibegin, iend, ifirst, iil, iindbl, iindw,
343 $ iindwk, iinfo, iinspl, iiu, ilast, in, indd,
344 $ inde2, inderr, indgp, indgrs, indwrk, itmp,
345 $ itmp2,
j, jblk, jj, liwmin, lwmin, nsplit,
346 $ nzcmin, offset, wbegin, wend
347 DOUBLE PRECISION bignum, cs, eps, pivmin, r1, r2, rmax, rmin,
348 $ rtol1, rtol2, safmin, scale, smlnum, sn,
349 $ thresh, tmp, tnrm, wl, wu
362 INTRINSIC max, min, sqrt
370 wantz =
lsame( jobz,
'V' )
371 alleig =
lsame( range,
'A' )
372 valeig =
lsame( range,
'V' )
373 indeig =
lsame( range,
'I' )
375 lquery = ( ( lwork.EQ.-1 ).OR.( liwork.EQ.-1 ) )
376 zquery = ( nzc.EQ.-1 )
402 ELSEIF( indeig )
THEN
409 IF( .NOT.( wantz .OR.
lsame( jobz,
'N' ) ) )
THEN
411 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) )
THEN
413 ELSE IF( n.LT.0 )
THEN
415 ELSE IF( valeig .AND. n.GT.0 .AND. wu.LE.wl )
THEN
417 ELSE IF( indeig .AND. ( iil.LT.1 .OR. iil.GT.n ) )
THEN
419 ELSE IF( indeig .AND. ( iiu.LT.iil .OR. iiu.GT.n ) )
THEN
421 ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
THEN
423 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery )
THEN
425 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery )
THEN
431 safmin =
dlamch(
'Safe minimum' )
432 eps =
dlamch(
'Precision' )
433 smlnum = safmin / eps
434 bignum = one / smlnum
435 rmin = sqrt( smlnum )
436 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
442 IF( wantz .AND. alleig )
THEN
444 ELSE IF( wantz .AND. valeig )
THEN
445 CALL
dlarrc(
'T', n, vl, vu, d, e, safmin,
446 $ nzcmin, itmp, itmp2, info )
447 ELSE IF( wantz .AND. indeig )
THEN
453 IF( zquery .AND. info.EQ.0 )
THEN
455 ELSE IF( nzc.LT.nzcmin .AND. .NOT.zquery )
THEN
462 CALL
xerbla(
'DSTEMR', -info )
465 ELSE IF( lquery .OR. zquery )
THEN
476 IF( alleig .OR. indeig )
THEN
480 IF( wl.LT.d( 1 ) .AND. wu.GE.d( 1 ) )
THEN
485 IF( wantz.AND.(.NOT.zquery) )
THEN
494 IF( .NOT.wantz )
THEN
495 CALL
dlae2( d(1), e(1), d(2), r1, r2 )
496 ELSE IF( wantz.AND.(.NOT.zquery) )
THEN
497 CALL
dlaev2( d(1), e(1), d(2), r1, r2, cs, sn )
500 $ (valeig.AND.(r2.GT.wl).AND.
502 $ (indeig.AND.(iil.EQ.1)) )
THEN
505 IF( wantz.AND.(.NOT.zquery) )
THEN
524 $ (valeig.AND.(r1.GT.wl).AND.
526 $ (indeig.AND.(iiu.EQ.2)) )
THEN
529 IF( wantz.AND.(.NOT.zquery) )
THEN
571 tnrm =
dlanst(
'M', n, d, e )
572 IF( tnrm.GT.zero .AND. tnrm.LT.rmin )
THEN
574 ELSE IF( tnrm.GT.rmax )
THEN
577 IF( scale.NE.one )
THEN
578 CALL
dscal( n, scale, d, 1 )
579 CALL
dscal( n-1, scale, e, 1 )
599 CALL
dlarrr( n, d, e, iinfo )
615 CALL
dcopy(n,d,1,work(indd),1)
619 work( inde2+
j-1 ) = e(
j)**2
623 IF( .NOT.wantz )
THEN
633 rtol2 = max( sqrt(eps)*5.0d-3, four * eps )
635 CALL
dlarre( range, n, wl, wu, iil, iiu, d, e,
636 $ work(inde2), rtol1, rtol2, thresh, nsplit,
637 $ iwork( iinspl ), m, w, work( inderr ),
638 $ work( indgp ), iwork( iindbl ),
639 $ iwork( iindw ), work( indgrs ), pivmin,
640 $ work( indwrk ), iwork( iindwk ), iinfo )
641 IF( iinfo.NE.0 )
THEN
642 info = 10 + abs( iinfo )
655 CALL
dlarrv( n, wl, wu, d, e,
656 $ pivmin, iwork( iinspl ), m,
657 $ 1, m, minrgp, rtol1, rtol2,
658 $ w, work( inderr ), work( indgp ), iwork( iindbl ),
659 $ iwork( iindw ), work( indgrs ), z, ldz,
660 $ isuppz, work( indwrk ), iwork( iindwk ), iinfo )
661 IF( iinfo.NE.0 )
THEN
662 info = 20 + abs( iinfo )
672 itmp = iwork( iindbl+
j-1 )
673 w(
j ) = w(
j ) + e( iwork( iinspl+itmp-1 ) )
683 DO 39 jblk = 1, iwork( iindbl+m-1 )
684 iend = iwork( iinspl+jblk-1 )
685 in = iend - ibegin + 1
690 IF( iwork( iindbl+wend ).EQ.jblk )
THEN
695 IF( wend.LT.wbegin )
THEN
700 offset = iwork(iindw+wbegin-1)-1
701 ifirst = iwork(iindw+wbegin-1)
702 ilast = iwork(iindw+wend-1)
705 $ work(indd+ibegin-1), work(inde2+ibegin-1),
706 $ ifirst, ilast, rtol2, offset, w(wbegin),
707 $ work( inderr+wbegin-1 ),
708 $ work( indwrk ), iwork( iindwk ), pivmin,
717 IF( scale.NE.one )
THEN
718 CALL
dscal( m, one / scale, w, 1 )
727 IF( nsplit.GT.1 .OR. n.EQ.2 )
THEN
728 IF( .NOT. wantz )
THEN
729 CALL
dlasrt(
'I', m, w, iinfo )
730 IF( iinfo.NE.0 )
THEN
739 IF( w( jj ).LT.tmp )
THEN
748 CALL
dswap( n, z( 1, i ), 1, z( 1,
j ), 1 )
749 itmp = isuppz( 2*i-1 )
750 isuppz( 2*i-1 ) = isuppz( 2*
j-1 )
751 isuppz( 2*
j-1 ) = itmp
753 isuppz( 2*i ) = isuppz( 2*
j )