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