1 SUBROUTINE ZLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK, 2 $ EPS3, SMLNUM, INFO ) 3* 4* -- LAPACK auxiliary routine (version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* September 30, 1994 8* 9* .. Scalar Arguments .. 10 LOGICAL NOINIT, RIGHTV 11 INTEGER INFO, LDB, LDH, N 12 DOUBLE PRECISION EPS3, SMLNUM 13 COMPLEX*16 W 14* .. 15* .. Array Arguments .. 16 DOUBLE PRECISION RWORK( * ) 17 COMPLEX*16 B( LDB, * ), H( LDH, * ), V( * ) 18* .. 19* 20* Purpose 21* ======= 22* 23* ZLAEIN uses inverse iteration to find a right or left eigenvector 24* corresponding to the eigenvalue W of a complex upper Hessenberg 25* matrix H. 26* 27* Arguments 28* ========= 29* 30* RIGHTV (input) LOGICAL 31* = .TRUE. : compute right eigenvector; 32* = .FALSE.: compute left eigenvector. 33* 34* NOINIT (input) LOGICAL 35* = .TRUE. : no initial vector supplied in V 36* = .FALSE.: initial vector supplied in V. 37* 38* N (input) INTEGER 39* The order of the matrix H. N >= 0. 40* 41* H (input) COMPLEX*16 array, dimension (LDH,N) 42* The upper Hessenberg matrix H. 43* 44* LDH (input) INTEGER 45* The leading dimension of the array H. LDH >= max(1,N). 46* 47* W (input) COMPLEX*16 48* The eigenvalue of H whose corresponding right or left 49* eigenvector is to be computed. 50* 51* V (input/output) COMPLEX*16 array, dimension (N) 52* On entry, if NOINIT = .FALSE., V must contain a starting 53* vector for inverse iteration; otherwise V need not be set. 54* On exit, V contains the computed eigenvector, normalized so 55* that the component of largest magnitude has magnitude 1; here 56* the magnitude of a complex number (x,y) is taken to be 57* |x| + |y|. 58* 59* B (workspace) COMPLEX*16 array, dimension (LDB,N) 60* 61* LDB (input) INTEGER 62* The leading dimension of the array B. LDB >= max(1,N). 63* 64* RWORK (workspace) DOUBLE PRECISION array, dimension (N) 65* 66* EPS3 (input) DOUBLE PRECISION 67* A small machine-dependent value which is used to perturb 68* close eigenvalues, and to replace zero pivots. 69* 70* SMLNUM (input) DOUBLE PRECISION 71* A machine-dependent value close to the underflow threshold. 72* 73* INFO (output) INTEGER 74* = 0: successful exit 75* = 1: inverse iteration did not converge; V is set to the 76* last iterate. 77* 78* ===================================================================== 79* 80* .. Parameters .. 81 DOUBLE PRECISION ONE, TENTH 82 PARAMETER ( ONE = 1.0D+0, TENTH = 1.0D-1 ) 83 COMPLEX*16 ZERO 84 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) 85* .. 86* .. Local Scalars .. 87 CHARACTER NORMIN, TRANS 88 INTEGER I, IERR, ITS, J 89 DOUBLE PRECISION GROWTO, NRMSML, ROOTN, RTEMP, SCALE, VNORM 90 COMPLEX*16 CDUM, EI, EJ, TEMP, X 91* .. 92* .. External Functions .. 93 INTEGER IZAMAX 94 DOUBLE PRECISION DZASUM, DZNRM2 95 COMPLEX*16 ZLADIV 96 EXTERNAL IZAMAX, DZASUM, DZNRM2, ZLADIV 97* .. 98* .. External Subroutines .. 99 EXTERNAL ZDSCAL, ZLATRS 100* .. 101* .. Intrinsic Functions .. 102 INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT 103* .. 104* .. Statement Functions .. 105 DOUBLE PRECISION CABS1 106* .. 107* .. Statement Function definitions .. 108 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) 109* .. 110* .. Executable Statements .. 111* 112 INFO = 0 113* 114* GROWTO is the threshold used in the acceptance test for an 115* eigenvector. 116* 117 ROOTN = SQRT( DBLE( N ) ) 118 GROWTO = TENTH / ROOTN 119 NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM 120* 121* Form B = H - W*I (except that the subdiagonal elements are not 122* stored). 123* 124 DO 20 J = 1, N 125 DO 10 I = 1, J - 1 126 B( I, J ) = H( I, J ) 127 10 CONTINUE 128 B( J, J ) = H( J, J ) - W 129 20 CONTINUE 130* 131 IF( NOINIT ) THEN 132* 133* Initialize V. 134* 135 DO 30 I = 1, N 136 V( I ) = EPS3 137 30 CONTINUE 138 ELSE 139* 140* Scale supplied initial vector. 141* 142 VNORM = DZNRM2( N, V, 1 ) 143 CALL ZDSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), V, 1 ) 144 END IF 145* 146 IF( RIGHTV ) THEN 147* 148* LU decomposition with partial pivoting of B, replacing zero 149* pivots by EPS3. 150* 151 DO 60 I = 1, N - 1 152 EI = H( I+1, I ) 153 IF( CABS1( B( I, I ) ).LT.CABS1( EI ) ) THEN 154* 155* Interchange rows and eliminate. 156* 157 X = ZLADIV( B( I, I ), EI ) 158 B( I, I ) = EI 159 DO 40 J = I + 1, N 160 TEMP = B( I+1, J ) 161 B( I+1, J ) = B( I, J ) - X*TEMP 162 B( I, J ) = TEMP 163 40 CONTINUE 164 ELSE 165* 166* Eliminate without interchange. 167* 168 IF( B( I, I ).EQ.ZERO ) 169 $ B( I, I ) = EPS3 170 X = ZLADIV( EI, B( I, I ) ) 171 IF( X.NE.ZERO ) THEN 172 DO 50 J = I + 1, N 173 B( I+1, J ) = B( I+1, J ) - X*B( I, J ) 174 50 CONTINUE 175 END IF 176 END IF 177 60 CONTINUE 178 IF( B( N, N ).EQ.ZERO ) 179 $ B( N, N ) = EPS3 180* 181 TRANS = 'N' 182* 183 ELSE 184* 185* UL decomposition with partial pivoting of B, replacing zero 186* pivots by EPS3. 187* 188 DO 90 J = N, 2, -1 189 EJ = H( J, J-1 ) 190 IF( CABS1( B( J, J ) ).LT.CABS1( EJ ) ) THEN 191* 192* Interchange columns and eliminate. 193* 194 X = ZLADIV( B( J, J ), EJ ) 195 B( J, J ) = EJ 196 DO 70 I = 1, J - 1 197 TEMP = B( I, J-1 ) 198 B( I, J-1 ) = B( I, J ) - X*TEMP 199 B( I, J ) = TEMP 200 70 CONTINUE 201 ELSE 202* 203* Eliminate without interchange. 204* 205 IF( B( J, J ).EQ.ZERO ) 206 $ B( J, J ) = EPS3 207 X = ZLADIV( EJ, B( J, J ) ) 208 IF( X.NE.ZERO ) THEN 209 DO 80 I = 1, J - 1 210 B( I, J-1 ) = B( I, J-1 ) - X*B( I, J ) 211 80 CONTINUE 212 END IF 213 END IF 214 90 CONTINUE 215 IF( B( 1, 1 ).EQ.ZERO ) 216 $ B( 1, 1 ) = EPS3 217* 218 TRANS = 'C' 219* 220 END IF 221* 222 NORMIN = 'N' 223 DO 110 ITS = 1, N 224* 225* Solve U*x = scale*v for a right eigenvector 226* or U'*x = scale*v for a left eigenvector, 227* overwriting x on v. 228* 229 CALL ZLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB, V, 230 $ SCALE, RWORK, IERR ) 231 NORMIN = 'Y' 232* 233* Test for sufficient growth in the norm of v. 234* 235 VNORM = DZASUM( N, V, 1 ) 236 IF( VNORM.GE.GROWTO*SCALE ) 237 $ GO TO 120 238* 239* Choose new orthogonal starting vector and try again. 240* 241 RTEMP = EPS3 / ( ROOTN+ONE ) 242 V( 1 ) = EPS3 243 DO 100 I = 2, N 244 V( I ) = RTEMP 245 100 CONTINUE 246 V( N-ITS+1 ) = V( N-ITS+1 ) - EPS3*ROOTN 247 110 CONTINUE 248* 249* Failure to find eigenvector in N iterations. 250* 251 INFO = 1 252* 253 120 CONTINUE 254* 255* Normalize eigenvector. 256* 257 I = IZAMAX( N, V, 1 ) 258 CALL ZDSCAL( N, ONE / CABS1( V( I ) ), V, 1 ) 259* 260 RETURN 261* 262* End of ZLAEIN 263* 264 END 265