1 *> \brief \b ZLAEIN 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 ZLAEIN + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaein.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaein.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaein.f">
21 * SUBROUTINE ZLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK,
22 * EPS3, SMLNUM, INFO )
24 * .. Scalar Arguments ..
25 * LOGICAL NOINIT, RIGHTV
26 * INTEGER INFO, LDB, LDH, N
27 * DOUBLE PRECISION EPS3, SMLNUM
30 * .. Array Arguments ..
31 * DOUBLE PRECISION RWORK( * )
32 * COMPLEX*16 B( LDB, * ), H( LDH, * ), V( * )
41 *> ZLAEIN uses inverse iteration to find a right or left eigenvector
42 *> corresponding to the eigenvalue W of a complex upper Hessenberg
52 *> = .TRUE. : compute right eigenvector;
53 *> = .FALSE.: compute left eigenvector.
59 *> = .TRUE. : no initial vector supplied in V
60 *> = .FALSE.: initial vector supplied in V.
66 *> The order of the matrix H. N >= 0.
71 *> H is COMPLEX*16 array, dimension (LDH,N)
72 *> The upper Hessenberg matrix H.
78 *> The leading dimension of the array H. LDH >= max(1,N).
84 *> The eigenvalue of H whose corresponding right or left
85 *> eigenvector is to be computed.
90 *> V is COMPLEX*16 array, dimension (N)
91 *> On entry, if NOINIT = .FALSE., V must contain a starting
92 *> vector for inverse iteration; otherwise V need not be set.
93 *> On exit, V contains the computed eigenvector, normalized so
94 *> that the component of largest magnitude has magnitude 1; here
95 *> the magnitude of a complex number (x,y) is taken to be
101 *> B is COMPLEX*16 array, dimension (LDB,N)
107 *> The leading dimension of the array B. LDB >= max(1,N).
112 *> RWORK is DOUBLE PRECISION array, dimension (N)
117 *> EPS3 is DOUBLE PRECISION
118 *> A small machine-dependent value which is used to perturb
119 *> close eigenvalues, and to replace zero pivots.
124 *> SMLNUM is DOUBLE PRECISION
125 *> A machine-dependent value close to the underflow threshold.
131 *> = 0: successful exit
132 *> = 1: inverse iteration did not converge; V is set to the
139 *> \author Univ. of Tennessee
140 *> \author Univ. of California Berkeley
141 *> \author Univ. of Colorado Denver
144 *> \date September 2012
146 *> \ingroup complex16OTHERauxiliary
148 * =====================================================================
149 SUBROUTINE ZLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK,
150 $ EPS3, SMLNUM, INFO )
152 * -- LAPACK auxiliary routine (version 3.4.2) --
153 * -- LAPACK is a software package provided by Univ. of Tennessee, --
154 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
157 * .. Scalar Arguments ..
158 LOGICAL NOINIT, RIGHTV
159 INTEGER INFO, LDB, LDH, N
160 DOUBLE PRECISION EPS3, SMLNUM
163 * .. Array Arguments ..
164 DOUBLE PRECISION RWORK( * )
165 COMPLEX*16 B( LDB, * ), H( LDH, * ), V( * )
168 * =====================================================================
171 DOUBLE PRECISION ONE, TENTH
172 PARAMETER ( ONE = 1.0D+0, TENTH = 1.0D-1 )
174 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
176 * .. Local Scalars ..
177 CHARACTER NORMIN, TRANS
178 INTEGER I, IERR, ITS, J
179 DOUBLE PRECISION GROWTO, NRMSML, ROOTN, RTEMP, SCALE, VNORM
180 COMPLEX*16 CDUM, EI, EJ, TEMP, X
182 * .. External Functions ..
184 DOUBLE PRECISION DZASUM, DZNRM2
186 EXTERNAL IZAMAX, DZASUM, DZNRM2, ZLADIV
188 * .. External Subroutines ..
189 EXTERNAL ZDSCAL, ZLATRS
191 * .. Intrinsic Functions ..
192 INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT
194 * .. Statement Functions ..
195 DOUBLE PRECISION CABS1
197 * .. Statement Function definitions ..
198 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
200 * .. Executable Statements ..
204 * GROWTO is the threshold used in the acceptance test for an
207 ROOTN = SQRT( DBLE( N ) )
208 GROWTO = TENTH / ROOTN
209 NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM
211 * Form B = H - W*I (except that the subdiagonal elements are not
216 B( I, J ) = H( I, J )
218 B( J, J ) = H( J, J ) - W
230 * Scale supplied initial vector.
232 VNORM = DZNRM2( N, V, 1 )
233 CALL ZDSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), V, 1 )
238 * LU decomposition with partial pivoting of B, replacing zero
243 IF( CABS1( B( I, I ) ).LT.CABS1( EI ) ) THEN
245 * Interchange rows and eliminate.
247 X = ZLADIV( B( I, I ), EI )
251 B( I+1, J ) = B( I, J ) - X*TEMP
256 * Eliminate without interchange.
258 IF( B( I, I ).EQ.ZERO )
260 X = ZLADIV( EI, B( I, I ) )
263 B( I+1, J ) = B( I+1, J ) - X*B( I, J )
268 IF( B( N, N ).EQ.ZERO )
275 * UL decomposition with partial pivoting of B, replacing zero
280 IF( CABS1( B( J, J ) ).LT.CABS1( EJ ) ) THEN
282 * Interchange columns and eliminate.
284 X = ZLADIV( B( J, J ), EJ )
288 B( I, J-1 ) = B( I, J ) - X*TEMP
293 * Eliminate without interchange.
295 IF( B( J, J ).EQ.ZERO )
297 X = ZLADIV( EJ, B( J, J ) )
300 B( I, J-1 ) = B( I, J-1 ) - X*B( I, J )
305 IF( B( 1, 1 ).EQ.ZERO )
315 * Solve U*x = scale*v for a right eigenvector
316 * or U**H *x = scale*v for a left eigenvector,
317 * overwriting x on v.
319 CALL ZLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB, V,
320 $ SCALE, RWORK, IERR )
323 * Test for sufficient growth in the norm of v.
325 VNORM = DZASUM( N, V, 1 )
326 IF( VNORM.GE.GROWTO*SCALE )
329 * Choose new orthogonal starting vector and try again.
331 RTEMP = EPS3 / ( ROOTN+ONE )
336 V( N-ITS+1 ) = V( N-ITS+1 ) - EPS3*ROOTN
339 * Failure to find eigenvector in N iterations.
345 * Normalize eigenvector.
347 I = IZAMAX( N, V, 1 )
348 CALL ZDSCAL( N, ONE / CABS1( V( I ) ), V, 1 )