172 SUBROUTINE dlaein( 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 DOUBLE PRECISION bignum, eps3, smlnum, wi, wr
186 DOUBLE PRECISION b( ldb, * ), h( ldh, * ), vi( * ), vr( * ),
193 DOUBLE PRECISION zero, one, tenth
194 parameter( zero = 0.0d+0, one = 1.0d+0, tenth = 1.0d-1 )
197 CHARACTER normin, trans
198 INTEGER i, i1, i2, i3, ierr, its,
j
199 DOUBLE PRECISION absbii, absbjj, ei, ej, growto, norm, nrmsml,
200 $ rec, rootn, scale, temp, vcrit, vmax, vnorm, w,
212 INTRINSIC abs, dble, max, sqrt
221 rootn = sqrt( dble( 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 =
dnrm2( n, vr, 1 )
251 CALL
dscal( 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
dlatrs(
'Upper', trans,
'Nonunit', normin, n,
b, ldb,
339 $ vr, scale, work, ierr )
344 vnorm =
dasum( 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
dscal( n, one / abs( vr( i ) ), vr, 1 )
385 rec = ( eps3*rootn ) / max( norm, nrmsml )
386 CALL
dscal( n, rec, vr, 1 )
387 CALL
dscal( n, rec, vi, 1 )
404 absbii =
dlapy2(
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 ) =
dasum( n-i,
b( i, i+1 ), ldb ) +
447 $
dasum( 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
dscal( n, rec, vr, 1 )
538 CALL
dscal( 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
dscal( n, rec, vr, 1 )
565 CALL
dscal( n, rec, vi, 1 )
575 CALL
dladiv( 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
dscal( n, one / vnorm, vr, 1 )
624 CALL
dscal( n, one / vnorm, vi, 1 )