c-----------------------------------------------------------------------
c\BeginDoc
c
c\Name: slaqrb
c
c\Description:
c Compute the eigenvalues and the Schur decomposition of an upper
c Hessenberg submatrix in rows and columns ILO to IHI. Only the
c last component of the Schur vectors are computed.
c
c This is mostly a modification of the LAPACK routine slahqr.
c
c\Usage:
c call slaqrb
c ( WANTT, N, ILO, IHI, H, LDH, WR, WI, Z, INFO )
c
c\Arguments
c WANTT Logical variable. (INPUT)
c = .TRUE. : the full Schur form T is required;
c = .FALSE.: only eigenvalues are required.
c
c N Integer. (INPUT)
c The order of the matrix H. N >= 0.
c
c ILO Integer. (INPUT)
c IHI Integer. (INPUT)
c It is assumed that H is already upper quasi-triangular in
c rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
c ILO = 1). SLAQRB works primarily with the Hessenberg
c submatrix in rows and columns ILO to IHI, but applies
c transformations to all of H if WANTT is .TRUE..
c 1 <= ILO <= max(1,IHI); IHI <= N.
c
c H Real array, dimension (LDH,N). (INPUT/OUTPUT)
c On entry, the upper Hessenberg matrix H.
c On exit, if WANTT is .TRUE., H is upper quasi-triangular in
c rows and columns ILO:IHI, with any 2-by-2 diagonal blocks in
c standard form. If WANTT is .FALSE., the contents of H are
c unspecified on exit.
c
c LDH Integer. (INPUT)
c The leading dimension of the array H. LDH >= max(1,N).
c
c WR Real array, dimension (N). (OUTPUT)
c WI Real array, dimension (N). (OUTPUT)
c The real and imaginary parts, respectively, of the computed
c eigenvalues ILO to IHI are stored in the corresponding
c elements of WR and WI. If two eigenvalues are computed as a
c complex conjugate pair, they are stored in consecutive
c elements of WR and WI, say the i-th and (i+1)th, with
c WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
c eigenvalues are stored in the same order as on the diagonal
c of the Schur form returned in H, with WR(i) = H(i,i), and, if
c H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
c WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
c
c Z Real array, dimension (N). (OUTPUT)
c On exit Z contains the last components of the Schur vectors.
c
c INFO Integer. (OUPUT)
c = 0: successful exit
c > 0: SLAQRB failed to compute all the eigenvalues ILO to IHI
c in a total of 30*(IHI-ILO+1) iterations; if INFO = i,
c elements i+1:ihi of WR and WI contain those eigenvalues
c which have been successfully computed.
c
c\Remarks
c 1. None.
c
c-----------------------------------------------------------------------
c
c\BeginLib
c
c\Local variables:
c xxxxxx real
c
c\Routines called:
c slabad LAPACK routine that computes machine constants.
c slamch LAPACK routine that determines machine constants.
c slanhs LAPACK routine that computes various norms of a matrix.
c slanv2 LAPACK routine that computes the Schur factorization of
c 2 by 2 nonsymmetric matrix in standard form.
c slarfg LAPACK Householder reflection construction routine.
c scopy Level 1 BLAS that copies one vector to another.
c srot Level 1 BLAS that applies a rotation to a 2 by 2 matrix.
c
c\Author
c Danny Sorensen Phuong Vu
c Richard Lehoucq CRPC / Rice University
c Dept. of Computational & Houston, Texas
c Applied Mathematics
c Rice University
c Houston, Texas
c
c\Revision history:
c xx/xx/92: Version ' 2.4'
c Modified from the LAPACK routine slahqr so that only the
c last component of the Schur vectors are computed.
c
c\SCCS Information: @(#)
c FILE: laqrb.F SID: 2.2 DATE OF SID: 8/27/96 RELEASE: 2
c
c\Remarks
c 1. None
c
c\EndLib
c
c-----------------------------------------------------------------------
c
subroutine slaqrb ( wantt, n, ilo, ihi, h, ldh, wr, wi,
& z, info )
c
c %------------------%
c | Scalar Arguments |
c %------------------%
c
logical wantt
integer ihi, ilo, info, ldh, n
c
c %-----------------%
c | Array Arguments |
c %-----------------%
c
Real
& h( ldh, * ), wi( * ), wr( * ), z( * )
c
c %------------%
c | Parameters |
c %------------%
c
Real
& zero, one, dat1, dat2
parameter (zero = 0.0E+0, one = 1.0E+0, dat1 = 7.5E-1,
& dat2 = -4.375E-1)
c
c %------------------------%
c | Local Scalars & Arrays |
c %------------------------%
c
integer i, i1, i2, itn, its, j, k, l, m, nh, nr
Real
& cs, h00, h10, h11, h12, h21, h22, h33, h33s,
& h43h34, h44, h44s, ovfl, s, smlnum, sn, sum,
& t1, t2, t3, tst1, ulp, unfl, v1, v2, v3
Real
& v( 3 ), work( 1 )
c
c %--------------------%
c | External Functions |
c %--------------------%
c
Real
& slamch, slanhs
external slamch, slanhs
c
c %----------------------%
c | External Subroutines |
c %----------------------%
c
external scopy, slabad, slanv2, slarfg, srot
c
c %-----------------------%
c | Executable Statements |
c %-----------------------%
c
info = 0
c
c %--------------------------%
c | Quick return if possible |
c %--------------------------%
c
if( n.eq.0 )
& return
if( ilo.eq.ihi ) then
wr( ilo ) = h( ilo, ilo )
wi( ilo ) = zero
return
end if
c
c %---------------------------------------------%
c | Initialize the vector of last components of |
c | the Schur vectors for accumulation. |
c %---------------------------------------------%
c
do 5 j = 1, n-1
z(j) = zero
5 continue
z(n) = one
c
nh = ihi - ilo + 1
c
c %-------------------------------------------------------------%
c | Set machine-dependent constants for the stopping criterion. |
c | If norm(H) <= sqrt(OVFL), overflow should not occur. |
c %-------------------------------------------------------------%
c
unfl = slamch( 'safe minimum' )
ovfl = one / unfl
call slabad( unfl, ovfl )
ulp = slamch( 'precision' )
smlnum = unfl*( nh / ulp )
c
c %---------------------------------------------------------------%
c | I1 and I2 are the indices of the first row and last column |
c | of H to which transformations must be applied. If eigenvalues |
c | only are computed, I1 and I2 are set inside the main loop. |
c | Zero out H(J+2,J) = ZERO for J=1:N if WANTT = .TRUE. |
c | else H(J+2,J) for J=ILO:IHI-ILO-1 if WANTT = .FALSE. |
c %---------------------------------------------------------------%
c
if( wantt ) then
i1 = 1
i2 = n
do 8 i=1,i2-2
h(i1+i+1,i) = zero
8 continue
else
do 9 i=1, ihi-ilo-1
h(ilo+i+1,ilo+i-1) = zero
9 continue
end if
c
c %---------------------------------------------------%
c | ITN is the total number of QR iterations allowed. |
c %---------------------------------------------------%
c
itn = 30*nh
c
c ------------------------------------------------------------------
c The main loop begins here. I is the loop index and decreases from
c IHI to ILO in steps of 1 or 2. Each iteration of the loop works
c with the active submatrix in rows and columns L to I.
c Eigenvalues I+1 to IHI have already converged. Either L = ILO or
c H(L,L-1) is negligible so that the matrix splits.
c ------------------------------------------------------------------
c
i = ihi
10 continue
l = ilo
if( i.lt.ilo )
& go to 150
c %--------------------------------------------------------------%
c | Perform QR iterations on rows and columns ILO to I until a |
c | submatrix of order 1 or 2 splits off at the bottom because a |
c | subdiagonal element has become negligible. |
c %--------------------------------------------------------------%
do 130 its = 0, itn
c
c %----------------------------------------------%
c | Look for a single small subdiagonal element. |
c %----------------------------------------------%
c
do 20 k = i, l + 1, -1
tst1 = abs( h( k-1, k-1 ) ) + abs( h( k, k ) )
if( tst1.eq.zero )
& tst1 = slanhs( '1', i-l+1, h( l, l ), ldh, work )
if( abs( h( k, k-1 ) ).le.max( ulp*tst1, smlnum ) )
& go to 30
20 continue
30 continue
l = k
if( l.gt.ilo ) then
c
c %------------------------%
c | H(L,L-1) is negligible |
c %------------------------%
c
h( l, l-1 ) = zero
end if
c
c %-------------------------------------------------------------%
c | Exit from loop if a submatrix of order 1 or 2 has split off |
c %-------------------------------------------------------------%
c
if( l.ge.i-1 )
& go to 140
c
c %---------------------------------------------------------%
c | Now the active submatrix is in rows and columns L to I. |
c | If eigenvalues only are being computed, only the active |
c | submatrix need be transformed. |
c %---------------------------------------------------------%
c
if( .not.wantt ) then
i1 = l
i2 = i
end if
c
if( its.eq.10 .or. its.eq.20 ) then
c
c %-------------------%
c | Exceptional shift |
c %-------------------%
c
s = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
h44 = dat1*s
h33 = h44
h43h34 = dat2*s*s
c
else
c
c %-----------------------------------------%
c | Prepare to use Wilkinson's double shift |
c %-----------------------------------------%
c
h44 = h( i, i )
h33 = h( i-1, i-1 )
h43h34 = h( i, i-1 )*h( i-1, i )
end if
c
c %-----------------------------------------------------%
c | Look for two consecutive small subdiagonal elements |
c %-----------------------------------------------------%
c
do 40 m = i - 2, l, -1
c
c %---------------------------------------------------------%
c | Determine the effect of starting the double-shift QR |
c | iteration at row M, and see if this would make H(M,M-1) |
c | negligible. |
c %---------------------------------------------------------%
c
h11 = h( m, m )
h22 = h( m+1, m+1 )
h21 = h( m+1, m )
h12 = h( m, m+1 )
h44s = h44 - h11
h33s = h33 - h11
v1 = ( h33s*h44s-h43h34 ) / h21 + h12
v2 = h22 - h11 - h33s - h44s
v3 = h( m+2, m+1 )
s = abs( v1 ) + abs( v2 ) + abs( v3 )
v1 = v1 / s
v2 = v2 / s
v3 = v3 / s
v( 1 ) = v1
v( 2 ) = v2
v( 3 ) = v3
if( m.eq.l )
& go to 50
h00 = h( m-1, m-1 )
h10 = h( m, m-1 )
tst1 = abs( v1 )*( abs( h00 )+abs( h11 )+abs( h22 ) )
if( abs( h10 )*( abs( v2 )+abs( v3 ) ).le.ulp*tst1 )
& go to 50
40 continue
50 continue
c
c %----------------------%
c | Double-shift QR step |
c %----------------------%
c
do 120 k = m, i - 1
c
c ------------------------------------------------------------
c The first iteration of this loop determines a reflection G
c from the vector V and applies it from left and right to H,
c thus creating a nonzero bulge below the subdiagonal.
c
c Each subsequent iteration determines a reflection G to
c restore the Hessenberg form in the (K-1)th column, and thus
c chases the bulge one step toward the bottom of the active
c submatrix. NR is the order of G.
c ------------------------------------------------------------
c
nr = min( 3, i-k+1 )
if( k.gt.m )
& call scopy( nr, h( k, k-1 ), 1, v, 1 )
call slarfg( nr, v( 1 ), v( 2 ), 1, t1 )
if( k.gt.m ) then
h( k, k-1 ) = v( 1 )
h( k+1, k-1 ) = zero
if( k.lt.i-1 )
& h( k+2, k-1 ) = zero
else if( m.gt.l ) then
h( k, k-1 ) = -h( k, k-1 )
end if
v2 = v( 2 )
t2 = t1*v2
if( nr.eq.3 ) then
v3 = v( 3 )
t3 = t1*v3
c
c %------------------------------------------------%
c | Apply G from the left to transform the rows of |
c | the matrix in columns K to I2. |
c %------------------------------------------------%
c
do 60 j = k, i2
sum = h( k, j ) + v2*h( k+1, j ) + v3*h( k+2, j )
h( k, j ) = h( k, j ) - sum*t1
h( k+1, j ) = h( k+1, j ) - sum*t2
h( k+2, j ) = h( k+2, j ) - sum*t3
60 continue
c
c %----------------------------------------------------%
c | Apply G from the right to transform the columns of |
c | the matrix in rows I1 to min(K+3,I). |
c %----------------------------------------------------%
c
do 70 j = i1, min( k+3, i )
sum = h( j, k ) + v2*h( j, k+1 ) + v3*h( j, k+2 )
h( j, k ) = h( j, k ) - sum*t1
h( j, k+1 ) = h( j, k+1 ) - sum*t2
h( j, k+2 ) = h( j, k+2 ) - sum*t3
70 continue
c
c %----------------------------------%
c | Accumulate transformations for Z |
c %----------------------------------%
c
sum = z( k ) + v2*z( k+1 ) + v3*z( k+2 )
z( k ) = z( k ) - sum*t1
z( k+1 ) = z( k+1 ) - sum*t2
z( k+2 ) = z( k+2 ) - sum*t3
else if( nr.eq.2 ) then
c
c %------------------------------------------------%
c | Apply G from the left to transform the rows of |
c | the matrix in columns K to I2. |
c %------------------------------------------------%
c
do 90 j = k, i2
sum = h( k, j ) + v2*h( k+1, j )
h( k, j ) = h( k, j ) - sum*t1
h( k+1, j ) = h( k+1, j ) - sum*t2
90 continue
c
c %----------------------------------------------------%
c | Apply G from the right to transform the columns of |
c | the matrix in rows I1 to min(K+3,I). |
c %----------------------------------------------------%
c
do 100 j = i1, i
sum = h( j, k ) + v2*h( j, k+1 )
h( j, k ) = h( j, k ) - sum*t1
h( j, k+1 ) = h( j, k+1 ) - sum*t2
100 continue
c
c %----------------------------------%
c | Accumulate transformations for Z |
c %----------------------------------%
c
sum = z( k ) + v2*z( k+1 )
z( k ) = z( k ) - sum*t1
z( k+1 ) = z( k+1 ) - sum*t2
end if
120 continue
130 continue
c
c %-------------------------------------------------------%
c | Failure to converge in remaining number of iterations |
c %-------------------------------------------------------%
c
info = i
return
140 continue
if( l.eq.i ) then
c
c %------------------------------------------------------%
c | H(I,I-1) is negligible: one eigenvalue has converged |
c %------------------------------------------------------%
c
wr( i ) = h( i, i )
wi( i ) = zero
else if( l.eq.i-1 ) then
c
c %--------------------------------------------------------%
c | H(I-1,I-2) is negligible; |
c | a pair of eigenvalues have converged. |
c | |
c | Transform the 2-by-2 submatrix to standard Schur form, |
c | and compute and store the eigenvalues. |
c %--------------------------------------------------------%
c
call slanv2( h( i-1, i-1 ), h( i-1, i ), h( i, i-1 ),
& h( i, i ), wr( i-1 ), wi( i-1 ), wr( i ), wi( i ),
& cs, sn )
if( wantt ) then
c
c %-----------------------------------------------------%
c | Apply the transformation to the rest of H and to Z, |
c | as required. |
c %-----------------------------------------------------%
c
if( i2.gt.i )
& call srot( i2-i, h( i-1, i+1 ), ldh, h( i, i+1 ), ldh,
& cs, sn )
call srot( i-i1-1, h( i1, i-1 ), 1, h( i1, i ), 1, cs, sn )
sum = cs*z( i-1 ) + sn*z( i )
z( i ) = cs*z( i ) - sn*z( i-1 )
z( i-1 ) = sum
end if
end if
c
c %---------------------------------------------------------%
c | Decrement number of remaining iterations, and return to |
c | start of the main loop with new value of I. |
c %---------------------------------------------------------%
c
itn = itn - its
i = l - 1
go to 10
150 continue
return
c
c %---------------%
c | End of slaqrb |
c %---------------%
c
end
