264 SUBROUTINE strsna( 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 REAL s( * ), sep( * ), t( ldt, * ), vl( ldvl, * ),
281 $ vr( ldvr, * ), work( ldwork, * )
288 parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
291 LOGICAL pair, somcon, wantbh, wants, wantsp
292 INTEGER i, ierr, ifst, ilst,
j, k, kase, ks, n2, nn
293 REAL bignum, cond, cs, delta, dumm, eps, est, lnrm,
294 $ mu, prod, prod1, prod2, rnrm, scale, smlnum, sn
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(
'STRSNA', -info )
383 IF( .NOT.
SELECT( 1 ) )
389 $ sep( 1 ) = abs( t( 1, 1 ) )
396 smlnum =
slamch(
'S' ) / eps
397 bignum = one / smlnum
398 CALL
slabad( 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 =
sdot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
439 rnrm =
snrm2( n, vr( 1, ks ), 1 )
440 lnrm =
snrm2( n, vl( 1, ks ), 1 )
441 s( ks ) = abs( prod ) / ( rnrm*lnrm )
446 prod1 =
sdot( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
447 prod1 = prod1 +
sdot( n, vr( 1, ks+1 ), 1, vl( 1, ks+1 ),
449 prod2 =
sdot( n, vl( 1, ks ), 1, vr( 1, ks+1 ), 1 )
450 prod2 = prod2 -
sdot( n, vl( 1, ks+1 ), 1, vr( 1, ks ),
453 $
snrm2( n, vr( 1, ks+1 ), 1 ) )
455 $
snrm2( n, vl( 1, ks+1 ), 1 ) )
456 cond =
slapy2( prod1, prod2 ) / ( rnrm*lnrm )
470 CALL
slacpy(
'Full', n, n, t, ldt, work, ldwork )
473 CALL
strexc(
'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 =
slapy2( 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
slacn2( nn, work( 1, n+2 ), work( 1, n+4 ), iwork,
549 CALL
slaqtr( .true., .true., n-1, work( 2, 2 ),
550 $ ldwork, dummy, dumm, scale,
551 $ work( 1, n+4 ), work( 1, n+6 ),
558 CALL
slaqtr( .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
slaqtr( .false., .true., n-1, work( 2, 2 ),
569 $ ldwork, dummy, dumm, scale,
570 $ work( 1, n+4 ), work( 1, n+6 ),
577 CALL
slaqtr( .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 strsna(JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK, INFO)
STRSNA
LOGICAL function lsame(CA, CB)
LSAME
REAL function slamch(CMACH)
SLAMCH
REAL function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine strexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, INFO)
STREXC
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine slacn2(N, V, X, ISGN, EST, KASE, ISAVE)
SLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine slaqtr(LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK, INFO)
SLAQTR solves a real quasi-triangular system of equations, or a complex quasi-triangular system of sp...
subroutine slabad(SMALL, LARGE)
SLABAD
REAL function sdot(N, SX, INCX, SY, INCY)
SDOT
REAL function snrm2(N, X, INCX)
SNRM2