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 )
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 ),
453 $
dnrm2( n, vr( 1, ks+1 ), 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 xerbla(SRNAME, INFO)
XERBLA
subroutine dlabad(SMALL, LARGE)
DLABAD
logical function lsame(CA, CB)
LSAME
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
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 dtrsna(JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK, INFO)
DTRSNA
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
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...
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
double precision function dnrm2(N, X, INCX)
DNRM2