172 SUBROUTINE slaein( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
173 $ ldb, work, eps3, smlnum, bignum, info )
181 LOGICAL noinit, rightv
182 INTEGER info, ldb, ldh, n
183 REAL bignum, eps3, smlnum, wi, wr
186 REAL b( ldb, * ), h( ldh, * ), vi( * ), vr( * ),
193 REAL zero, one, tenth
194 parameter( zero = 0.0e+0, one = 1.0e+0, tenth = 1.0e-1 )
197 CHARACTER normin, trans
198 INTEGER i, i1, i2, i3, ierr, its,
j
199 REAL absbii, absbjj, ei, ej, growto, norm, nrmsml,
200 $ rec, rootn, scale, temp, vcrit, vmax, vnorm, w,
212 INTRINSIC abs, max,
REAL, sqrt
221 rootn = sqrt(
REAL( N ) )
222 growto = tenth / rootn
223 nrmsml = max( one, eps3*rootn )*smlnum
230 b( i,
j ) = h( i,
j )
232 b(
j,
j ) = h(
j,
j ) - wr
235 IF( wi.EQ.zero )
THEN
250 vnorm =
snrm2( n, vr, 1 )
251 CALL
sscal( n, ( eps3*rootn ) / max( vnorm, nrmsml ), vr,
262 IF( abs(
b( i, i ) ).LT.abs( ei ) )
THEN
270 b( i+1,
j ) =
b( i,
j ) -
x*temp
277 IF(
b( i, i ).EQ.zero )
282 b( i+1,
j ) =
b( i+1,
j ) -
x*
b( i,
j )
287 IF(
b( n, n ).EQ.zero )
299 IF( abs(
b(
j,
j ) ).LT.abs( ej ) )
THEN
307 b( i,
j-1 ) =
b( i,
j ) -
x*temp
314 IF(
b(
j,
j ).EQ.zero )
319 b( i,
j-1 ) =
b( i,
j-1 ) -
x*
b( i,
j )
324 IF(
b( 1, 1 ).EQ.zero )
338 CALL
slatrs(
'Upper', trans,
'Nonunit', normin, n,
b, ldb,
339 $ vr, scale, work, ierr )
344 vnorm =
sasum( n, vr, 1 )
345 IF( vnorm.GE.growto*scale )
350 temp = eps3 / ( rootn+one )
355 vr( n-its+1 ) = vr( n-its+1 ) - eps3*rootn
367 CALL
sscal( n, one / abs( vr( i ) ), vr, 1 )
385 rec = ( eps3*rootn ) / max( norm, nrmsml )
386 CALL
sscal( n, rec, vr, 1 )
387 CALL
sscal( n, rec, vi, 1 )
404 absbii =
slapy2(
b( i, i ),
b( i+1, i ) )
406 IF( absbii.LT.abs( ei ) )
THEN
411 xi =
b( i+1, i ) / ei
416 b( i+1,
j ) =
b( i,
j ) - xr*temp
417 b(
j+1, i+1 ) =
b(
j+1, i ) - xi*temp
422 b( i+1, i+1 ) =
b( i+1, i+1 ) - xi*wi
423 b( i+2, i+1 ) =
b( i+2, i+1 ) + xr*wi
428 IF( absbii.EQ.zero )
THEN
433 ei = ( ei / absbii ) / absbii
437 b( i+1,
j ) =
b( i+1,
j ) - xr*
b( i,
j ) +
439 b(
j+1, i+1 ) = -xr*
b(
j+1, i ) - xi*
b( i,
j )
441 b( i+2, i+1 ) =
b( i+2, i+1 ) - wi
446 work( i ) =
sasum( n-i,
b( i, i+1 ), ldb ) +
447 $
sasum( n-i,
b( i+2, i ), 1 )
449 IF(
b( n, n ).EQ.zero .AND.
b( n+1, n ).EQ.zero )
472 IF( absbjj.LT.abs( ej ) )
THEN
477 xi =
b(
j+1,
j ) / ej
482 b( i,
j-1 ) =
b( i,
j ) - xr*temp
483 b(
j, i ) =
b(
j+1, i ) - xi*temp
488 b(
j-1,
j-1 ) =
b(
j-1,
j-1 ) + xi*wi
489 b(
j,
j-1 ) =
b(
j,
j-1 ) - xr*wi
494 IF( absbjj.EQ.zero )
THEN
499 ej = ( ej / absbjj ) / absbjj
503 b( i,
j-1 ) =
b( i,
j-1 ) - xr*
b( i,
j ) +
505 b(
j, i ) = -xr*
b(
j+1, i ) - xi*
b( i,
j )
507 b(
j,
j-1 ) =
b(
j,
j-1 ) + wi
515 IF(
b( 1, 1 ).EQ.zero .AND.
b( 2, 1 ).EQ.zero )
533 DO 250 i = i1, i2, i3
535 IF( work( i ).GT.vcrit )
THEN
537 CALL
sscal( n, rec, vr, 1 )
538 CALL
sscal( n, rec, vi, 1 )
548 xr = xr -
b( i,
j )*vr(
j ) +
b(
j+1, i )*vi(
j )
549 xi = xi -
b( i,
j )*vi(
j ) -
b(
j+1, i )*vr(
j )
553 xr = xr -
b(
j, i )*vr(
j ) +
b( i+1,
j )*vi(
j )
554 xi = xi -
b(
j, i )*vi(
j ) -
b( i+1,
j )*vr(
j )
558 w = abs(
b( i, i ) ) + abs(
b( i+1, i ) )
559 IF( w.GT.smlnum )
THEN
561 w1 = abs( xr ) + abs( xi )
562 IF( w1.GT.w*bignum )
THEN
564 CALL
sscal( n, rec, vr, 1 )
565 CALL
sscal( n, rec, vi, 1 )
575 CALL
sladiv( xr, xi,
b( i, i ),
b( i+1, i ), vr( i ),
577 vmax = max( abs( vr( i ) )+abs( vi( i ) ), vmax )
578 vcrit = bignum / vmax
595 IF( vnorm.GE.growto*scale )
600 y = eps3 / ( rootn+one )
608 vr( n-its+1 ) = vr( n-its+1 ) - eps3*rootn
621 vnorm = max( vnorm, abs( vr( i ) )+abs( vi( i ) ) )
623 CALL
sscal( n, one / vnorm, vr, 1 )
624 CALL
sscal( n, one / vnorm, vi, 1 )
subroutine sladiv(A, B, C, D, P, Q)
SLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
INTEGER function isamax(N, SX, INCX)
ISAMAX
REAL function slapy2(X, Y)
SLAPY2 returns sqrt(x2+y2).
subroutine slaein(RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B, LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO)
SLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iterat...
subroutine slatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
SLATRS solves a triangular system of equations with the scale factor set to prevent overflow...
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
REAL function sasum(N, SX, INCX)
SASUM
REAL function snrm2(N, X, INCX)
SNRM2
subroutine sscal(N, SA, SX, INCX)
SSCAL