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 )
subroutine dlaein(RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B, LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO)
DLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iterat...
subroutine dlatrs(UPLO, TRANS, DIAG, NORMIN, N, A, LDA, X, SCALE, CNORM, INFO)
DLATRS 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
subroutine dscal(N, DA, DX, INCX)
DSCAL
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine dladiv(A, B, C, D, P, Q)
DLADIV performs complex division in real arithmetic, avoiding unnecessary overflow.
integer function idamax(N, DX, INCX)
IDAMAX
double precision function dnrm2(N, X, INCX)
DNRM2
double precision function dasum(N, DX, INCX)
DASUM