264 SUBROUTINE dtrsna( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
265 $ ldvr, s, sep, mm, m, work, ldwork, iwork,
274 CHARACTER HOWMNY, JOB
275 INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
280 DOUBLE PRECISION S( * ), SEP( * ), T( ldt, * ), VL( ldvl, * ),
281 $ vr( ldvr, * ), work( ldwork, * )
287 DOUBLE PRECISION ZERO, ONE, TWO
288 parameter ( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
291 LOGICAL PAIR, SOMCON, WANTBH, WANTS, WANTSP
292 INTEGER I, IERR, IFST, ILST, J, K, KASE, KS, N2, NN
293 DOUBLE PRECISION BIGNUM, COND, CS, DELTA, DUMM, EPS, EST, LNRM,
294 $ mu, prod, prod1, prod2, rnrm, scale, smlnum, sn
298 DOUBLE PRECISION DUMMY( 1 )
302 DOUBLE PRECISION DDOT, DLAMCH, DLAPY2, DNRM2
303 EXTERNAL lsame, ddot, dlamch, dlapy2, dnrm2
309 INTRINSIC abs, max, sqrt
315 wantbh = lsame( job,
'B' )
316 wants = lsame( job,
'E' ) .OR. wantbh
317 wantsp = lsame( job,
'V' ) .OR. wantbh
319 somcon = lsame( howmny,
'S' )
322 IF( .NOT.wants .AND. .NOT.wantsp )
THEN
324 ELSE IF( .NOT.lsame( howmny,
'A' ) .AND. .NOT.somcon )
THEN
326 ELSE IF( n.LT.0 )
THEN
328 ELSE IF( ldt.LT.max( 1, n ) )
THEN
330 ELSE IF( ldvl.LT.1 .OR. ( wants .AND. ldvl.LT.n ) )
THEN
332 ELSE IF( ldvr.LT.1 .OR. ( wants .AND. ldvr.LT.n ) )
THEN
347 IF( t( k+1, k ).EQ.zero )
THEN
352 IF(
SELECT( k ) .OR.
SELECT( k+1 ) )
367 ELSE IF( ldwork.LT.1 .OR. ( wantsp .AND. ldwork.LT.n ) )
THEN
372 CALL xerbla(
'DTRSNA', -info )
383 IF( .NOT.
SELECT( 1 ) )
389 $ sep( 1 ) = abs( t( 1, 1 ) )
396 smlnum = dlamch(
'S' ) / eps
397 bignum = one / smlnum
398 CALL dlabad( smlnum, bignum )
411 $ pair = t( k+1, k ).NE.zero
419 IF( .NOT.
SELECT( k ) .AND. .NOT.
SELECT( k+1 ) )
422 IF( .NOT.
SELECT( k ) )
438 prod = ddot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
439 rnrm = dnrm2( n, vr( 1, ks ), 1 )
440 lnrm = dnrm2( n, vl( 1, ks ), 1 )
441 s( ks ) = abs( prod ) / ( rnrm*lnrm )
446 prod1 = ddot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
447 prod1 = prod1 + ddot( n, vr( 1, ks+1 ), 1, vl( 1, ks+1 ),
449 prod2 = ddot( n, vl( 1, ks ), 1, vr( 1, ks+1 ), 1 )
450 prod2 = prod2 - ddot( n, vl( 1, ks+1 ), 1, vr( 1, ks ),
452 rnrm = dlapy2( dnrm2( n, vr( 1, ks ), 1 ),
453 $ dnrm2( n, vr( 1, ks+1 ), 1 ) )
454 lnrm = dlapy2( dnrm2( n, vl( 1, ks ), 1 ),
455 $ dnrm2( n, vl( 1, ks+1 ), 1 ) )
456 cond = dlapy2( prod1, prod2 ) / ( rnrm*lnrm )
470 CALL dlacpy(
'Full', n, n, t, ldt, work, ldwork )
473 CALL dtrexc(
'No Q', n, work, ldwork, dummy, 1, ifst, ilst,
474 $ work( 1, n+1 ), ierr )
476 IF( ierr.EQ.1 .OR. ierr.EQ.2 )
THEN
486 IF( work( 2, 1 ).EQ.zero )
THEN
491 work( i, i ) = work( i, i ) - work( 1, 1 )
505 mu = sqrt( abs( work( 1, 2 ) ) )*
506 $ sqrt( abs( work( 2, 1 ) ) )
507 delta = dlapy2( mu, work( 2, 1 ) )
509 sn = -work( 2, 1 ) / delta
523 work( 2, j ) = cs*work( 2, j )
524 work( j, j ) = work( j, j ) - work( 1, 1 )
528 work( 1, n+1 ) = two*mu
530 work( i, n+1 ) = sn*work( 1, i+1 )
541 CALL dlacn2( nn, work( 1, n+2 ), work( 1, n+4 ), iwork,
549 CALL dlaqtr( .true., .true., n-1, work( 2, 2 ),
550 $ ldwork, dummy, dumm, scale,
551 $ work( 1, n+4 ), work( 1, n+6 ),
558 CALL dlaqtr( .true., .false., n-1, work( 2, 2 ),
559 $ ldwork, work( 1, n+1 ), mu, scale,
560 $ work( 1, n+4 ), work( 1, n+6 ),
568 CALL dlaqtr( .false., .true., n-1, work( 2, 2 ),
569 $ ldwork, dummy, dumm, scale,
570 $ work( 1, n+4 ), work( 1, n+6 ),
577 CALL dlaqtr( .false., .false., n-1,
578 $ work( 2, 2 ), ldwork,
579 $ work( 1, n+1 ), mu, scale,
580 $ work( 1, n+4 ), work( 1, n+6 ),
590 sep( ks ) = scale / max( est, smlnum )
592 $ sep( ks+1 ) = sep( ks )
subroutine dtrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, INFO)
DTREXC
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dtrsna(JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK, INFO)
DTRSNA
subroutine dlaqtr(LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK, INFO)
DLAQTR solves a real quasi-triangular system of equations, or a complex quasi-triangular system of sp...
subroutine dlacn2(N, V, X, ISGN, EST, KASE, ISAVE)
DLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...