1 *> \brief \b SLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration.
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download SLAEIN + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaein.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaein.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaein.f">
21 * SUBROUTINE SLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
22 * LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO )
24 * .. Scalar Arguments ..
25 * LOGICAL NOINIT, RIGHTV
26 * INTEGER INFO, LDB, LDH, N
27 * REAL BIGNUM, EPS3, SMLNUM, WI, WR
29 * .. Array Arguments ..
30 * REAL B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
40 *> SLAEIN uses inverse iteration to find a right or left eigenvector
41 *> corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg
51 *> = .TRUE. : compute right eigenvector;
52 *> = .FALSE.: compute left eigenvector.
58 *> = .TRUE. : no initial vector supplied in (VR,VI).
59 *> = .FALSE.: initial vector supplied in (VR,VI).
65 *> The order of the matrix H. N >= 0.
70 *> H is REAL array, dimension (LDH,N)
71 *> The upper Hessenberg matrix H.
77 *> The leading dimension of the array H. LDH >= max(1,N).
88 *> The real and imaginary parts of the eigenvalue of H whose
89 *> corresponding right or left eigenvector is to be computed.
94 *> VR is REAL array, dimension (N)
99 *> VI is REAL array, dimension (N)
100 *> On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain
101 *> a real starting vector for inverse iteration using the real
102 *> eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI
103 *> must contain the real and imaginary parts of a complex
104 *> starting vector for inverse iteration using the complex
105 *> eigenvalue (WR,WI); otherwise VR and VI need not be set.
106 *> On exit, if WI = 0.0 (real eigenvalue), VR contains the
107 *> computed real eigenvector; if WI.ne.0.0 (complex eigenvalue),
108 *> VR and VI contain the real and imaginary parts of the
109 *> computed complex eigenvector. The eigenvector is normalized
110 *> so that the component of largest magnitude has magnitude 1;
111 *> here the magnitude of a complex number (x,y) is taken to be
113 *> VI is not referenced if WI = 0.0.
118 *> B is REAL array, dimension (LDB,N)
124 *> The leading dimension of the array B. LDB >= N+1.
129 *> WORK is REAL array, dimension (N)
135 *> A small machine-dependent value which is used to perturb
136 *> close eigenvalues, and to replace zero pivots.
142 *> A machine-dependent value close to the underflow threshold.
148 *> A machine-dependent value close to the overflow threshold.
154 *> = 0: successful exit
155 *> = 1: inverse iteration did not converge; VR is set to the
156 *> last iterate, and so is VI if WI.ne.0.0.
162 *> \author Univ. of Tennessee
163 *> \author Univ. of California Berkeley
164 *> \author Univ. of Colorado Denver
167 *> \date September 2012
169 *> \ingroup realOTHERauxiliary
171 * =====================================================================
172 SUBROUTINE SLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
173 $ LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO )
175 * -- LAPACK auxiliary routine (version 3.4.2) --
176 * -- LAPACK is a software package provided by Univ. of Tennessee, --
177 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
180 * .. Scalar Arguments ..
181 LOGICAL NOINIT, RIGHTV
182 INTEGER INFO, LDB, LDH, N
183 REAL BIGNUM, EPS3, SMLNUM, WI, WR
185 * .. Array Arguments ..
186 REAL B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
190 * =====================================================================
193 REAL ZERO, ONE, TENTH
194 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TENTH = 1.0E-1 )
196 * .. Local Scalars ..
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,
203 * .. External Functions ..
205 REAL SASUM, SLAPY2, SNRM2
206 EXTERNAL ISAMAX, SASUM, SLAPY2, SNRM2
208 * .. External Subroutines ..
209 EXTERNAL SLADIV, SLATRS, SSCAL
211 * .. Intrinsic Functions ..
212 INTRINSIC ABS, MAX, REAL, SQRT
214 * .. Executable Statements ..
218 * GROWTO is the threshold used in the acceptance test for an
221 ROOTN = SQRT( REAL( N ) )
222 GROWTO = TENTH / ROOTN
223 NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM
225 * Form B = H - (WR,WI)*I (except that the subdiagonal elements and
226 * the imaginary parts of the diagonal elements are not stored).
230 B( I, J ) = H( I, J )
232 B( J, J ) = H( J, J ) - WR
235 IF( WI.EQ.ZERO ) THEN
241 * Set initial vector.
248 * Scale supplied initial vector.
250 VNORM = SNRM2( N, VR, 1 )
251 CALL SSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), VR,
257 * LU decomposition with partial pivoting of B, replacing zero
262 IF( ABS( B( I, I ) ).LT.ABS( EI ) ) THEN
264 * Interchange rows and eliminate.
270 B( I+1, J ) = B( I, J ) - X*TEMP
275 * Eliminate without interchange.
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 )
294 * UL decomposition with partial pivoting of B, replacing zero
299 IF( ABS( B( J, J ) ).LT.ABS( EJ ) ) THEN
301 * Interchange columns and eliminate.
307 B( I, J-1 ) = B( I, J ) - X*TEMP
312 * Eliminate without interchange.
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 )
334 * Solve U*x = scale*v for a right eigenvector
335 * or U**T*x = scale*v for a left eigenvector,
336 * overwriting x on v.
338 CALL SLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB,
339 $ VR, SCALE, WORK, IERR )
342 * Test for sufficient growth in the norm of v.
344 VNORM = SASUM( N, VR, 1 )
345 IF( VNORM.GE.GROWTO*SCALE )
348 * Choose new orthogonal starting vector and try again.
350 TEMP = EPS3 / ( ROOTN+ONE )
355 VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
358 * Failure to find eigenvector in N iterations.
364 * Normalize eigenvector.
366 I = ISAMAX( N, VR, 1 )
367 CALL SSCAL( N, ONE / ABS( VR( I ) ), VR, 1 )
370 * Complex eigenvalue.
374 * Set initial vector.
382 * Scale supplied initial vector.
384 NORM = SLAPY2( SNRM2( N, VR, 1 ), SNRM2( N, VI, 1 ) )
385 REC = ( EPS3*ROOTN ) / MAX( NORM, NRMSML )
386 CALL SSCAL( N, REC, VR, 1 )
387 CALL SSCAL( N, REC, VI, 1 )
392 * LU decomposition with partial pivoting of B, replacing zero
395 * The imaginary part of the (i,j)-th element of U is stored in
404 ABSBII = SLAPY2( B( I, I ), B( I+1, I ) )
406 IF( ABSBII.LT.ABS( EI ) ) THEN
408 * Interchange rows and eliminate.
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
426 * Eliminate without interchanging rows.
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
444 * Compute 1-norm of offdiagonal elements of i-th row.
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 )
458 * UL decomposition with partial pivoting of conjg(B),
459 * replacing zero pivots by EPS3.
461 * The imaginary part of the (i,j)-th element of U is stored in
471 ABSBJJ = SLAPY2( B( J, J ), B( J+1, J ) )
472 IF( ABSBJJ.LT.ABS( EJ ) ) THEN
474 * Interchange columns and eliminate
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
492 * Eliminate without interchange.
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
510 * Compute 1-norm of offdiagonal elements of j-th column.
512 WORK( J ) = SASUM( J-1, B( 1, J ), 1 ) +
513 $ SASUM( J-1, B( J+1, 1 ), LDB )
515 IF( B( 1, 1 ).EQ.ZERO .AND. B( 2, 1 ).EQ.ZERO )
529 * Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector,
530 * or U**T*(xr,xi) = scale*(vr,vi) for a left eigenvector,
531 * overwriting (xr,xi) on (vr,vi).
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 )
573 * Divide by diagonal element of B.
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
592 * Test for sufficient growth in the norm of (VR,VI).
594 VNORM = SASUM( N, VR, 1 ) + SASUM( N, VI, 1 )
595 IF( VNORM.GE.GROWTO*SCALE )
598 * Choose a new orthogonal starting vector and try again.
600 Y = EPS3 / ( ROOTN+ONE )
608 VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
611 * Failure to find eigenvector in N iterations
617 * Normalize eigenvector.
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 )