1*> \brief \b CLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download CLAEIN + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claein.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claein.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claein.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE CLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK, 22* EPS3, SMLNUM, INFO ) 23* 24* .. Scalar Arguments .. 25* LOGICAL NOINIT, RIGHTV 26* INTEGER INFO, LDB, LDH, N 27* REAL EPS3, SMLNUM 28* COMPLEX W 29* .. 30* .. Array Arguments .. 31* REAL RWORK( * ) 32* COMPLEX B( LDB, * ), H( LDH, * ), V( * ) 33* .. 34* 35* 36*> \par Purpose: 37* ============= 38*> 39*> \verbatim 40*> 41*> CLAEIN uses inverse iteration to find a right or left eigenvector 42*> corresponding to the eigenvalue W of a complex upper Hessenberg 43*> matrix H. 44*> \endverbatim 45* 46* Arguments: 47* ========== 48* 49*> \param[in] RIGHTV 50*> \verbatim 51*> RIGHTV is LOGICAL 52*> = .TRUE. : compute right eigenvector; 53*> = .FALSE.: compute left eigenvector. 54*> \endverbatim 55*> 56*> \param[in] NOINIT 57*> \verbatim 58*> NOINIT is LOGICAL 59*> = .TRUE. : no initial vector supplied in V 60*> = .FALSE.: initial vector supplied in V. 61*> \endverbatim 62*> 63*> \param[in] N 64*> \verbatim 65*> N is INTEGER 66*> The order of the matrix H. N >= 0. 67*> \endverbatim 68*> 69*> \param[in] H 70*> \verbatim 71*> H is COMPLEX array, dimension (LDH,N) 72*> The upper Hessenberg matrix H. 73*> \endverbatim 74*> 75*> \param[in] LDH 76*> \verbatim 77*> LDH is INTEGER 78*> The leading dimension of the array H. LDH >= max(1,N). 79*> \endverbatim 80*> 81*> \param[in] W 82*> \verbatim 83*> W is COMPLEX 84*> The eigenvalue of H whose corresponding right or left 85*> eigenvector is to be computed. 86*> \endverbatim 87*> 88*> \param[in,out] V 89*> \verbatim 90*> V is COMPLEX 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 96*> |x| + |y|. 97*> \endverbatim 98*> 99*> \param[out] B 100*> \verbatim 101*> B is COMPLEX array, dimension (LDB,N) 102*> \endverbatim 103*> 104*> \param[in] LDB 105*> \verbatim 106*> LDB is INTEGER 107*> The leading dimension of the array B. LDB >= max(1,N). 108*> \endverbatim 109*> 110*> \param[out] RWORK 111*> \verbatim 112*> RWORK is REAL array, dimension (N) 113*> \endverbatim 114*> 115*> \param[in] EPS3 116*> \verbatim 117*> EPS3 is REAL 118*> A small machine-dependent value which is used to perturb 119*> close eigenvalues, and to replace zero pivots. 120*> \endverbatim 121*> 122*> \param[in] SMLNUM 123*> \verbatim 124*> SMLNUM is REAL 125*> A machine-dependent value close to the underflow threshold. 126*> \endverbatim 127*> 128*> \param[out] INFO 129*> \verbatim 130*> INFO is INTEGER 131*> = 0: successful exit 132*> = 1: inverse iteration did not converge; V is set to the 133*> last iterate. 134*> \endverbatim 135* 136* Authors: 137* ======== 138* 139*> \author Univ. of Tennessee 140*> \author Univ. of California Berkeley 141*> \author Univ. of Colorado Denver 142*> \author NAG Ltd. 143* 144*> \ingroup complexOTHERauxiliary 145* 146* ===================================================================== 147 SUBROUTINE CLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK, 148 $ EPS3, SMLNUM, INFO ) 149* 150* -- LAPACK auxiliary routine -- 151* -- LAPACK is a software package provided by Univ. of Tennessee, -- 152* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 153* 154* .. Scalar Arguments .. 155 LOGICAL NOINIT, RIGHTV 156 INTEGER INFO, LDB, LDH, N 157 REAL EPS3, SMLNUM 158 COMPLEX W 159* .. 160* .. Array Arguments .. 161 REAL RWORK( * ) 162 COMPLEX B( LDB, * ), H( LDH, * ), V( * ) 163* .. 164* 165* ===================================================================== 166* 167* .. Parameters .. 168 REAL ONE, TENTH 169 PARAMETER ( ONE = 1.0E+0, TENTH = 1.0E-1 ) 170 COMPLEX ZERO 171 PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ) ) 172* .. 173* .. Local Scalars .. 174 CHARACTER NORMIN, TRANS 175 INTEGER I, IERR, ITS, J 176 REAL GROWTO, NRMSML, ROOTN, RTEMP, SCALE, VNORM 177 COMPLEX CDUM, EI, EJ, TEMP, X 178* .. 179* .. External Functions .. 180 INTEGER ICAMAX 181 REAL SCASUM, SCNRM2 182 COMPLEX CLADIV 183 EXTERNAL ICAMAX, SCASUM, SCNRM2, CLADIV 184* .. 185* .. External Subroutines .. 186 EXTERNAL CLATRS, CSSCAL 187* .. 188* .. Intrinsic Functions .. 189 INTRINSIC ABS, AIMAG, MAX, REAL, SQRT 190* .. 191* .. Statement Functions .. 192 REAL CABS1 193* .. 194* .. Statement Function definitions .. 195 CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) ) 196* .. 197* .. Executable Statements .. 198* 199 INFO = 0 200* 201* GROWTO is the threshold used in the acceptance test for an 202* eigenvector. 203* 204 ROOTN = SQRT( REAL( N ) ) 205 GROWTO = TENTH / ROOTN 206 NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM 207* 208* Form B = H - W*I (except that the subdiagonal elements are not 209* stored). 210* 211 DO 20 J = 1, N 212 DO 10 I = 1, J - 1 213 B( I, J ) = H( I, J ) 214 10 CONTINUE 215 B( J, J ) = H( J, J ) - W 216 20 CONTINUE 217* 218 IF( NOINIT ) THEN 219* 220* Initialize V. 221* 222 DO 30 I = 1, N 223 V( I ) = EPS3 224 30 CONTINUE 225 ELSE 226* 227* Scale supplied initial vector. 228* 229 VNORM = SCNRM2( N, V, 1 ) 230 CALL CSSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), V, 1 ) 231 END IF 232* 233 IF( RIGHTV ) THEN 234* 235* LU decomposition with partial pivoting of B, replacing zero 236* pivots by EPS3. 237* 238 DO 60 I = 1, N - 1 239 EI = H( I+1, I ) 240 IF( CABS1( B( I, I ) ).LT.CABS1( EI ) ) THEN 241* 242* Interchange rows and eliminate. 243* 244 X = CLADIV( B( I, I ), EI ) 245 B( I, I ) = EI 246 DO 40 J = I + 1, N 247 TEMP = B( I+1, J ) 248 B( I+1, J ) = B( I, J ) - X*TEMP 249 B( I, J ) = TEMP 250 40 CONTINUE 251 ELSE 252* 253* Eliminate without interchange. 254* 255 IF( B( I, I ).EQ.ZERO ) 256 $ B( I, I ) = EPS3 257 X = CLADIV( EI, B( I, I ) ) 258 IF( X.NE.ZERO ) THEN 259 DO 50 J = I + 1, N 260 B( I+1, J ) = B( I+1, J ) - X*B( I, J ) 261 50 CONTINUE 262 END IF 263 END IF 264 60 CONTINUE 265 IF( B( N, N ).EQ.ZERO ) 266 $ B( N, N ) = EPS3 267* 268 TRANS = 'N' 269* 270 ELSE 271* 272* UL decomposition with partial pivoting of B, replacing zero 273* pivots by EPS3. 274* 275 DO 90 J = N, 2, -1 276 EJ = H( J, J-1 ) 277 IF( CABS1( B( J, J ) ).LT.CABS1( EJ ) ) THEN 278* 279* Interchange columns and eliminate. 280* 281 X = CLADIV( B( J, J ), EJ ) 282 B( J, J ) = EJ 283 DO 70 I = 1, J - 1 284 TEMP = B( I, J-1 ) 285 B( I, J-1 ) = B( I, J ) - X*TEMP 286 B( I, J ) = TEMP 287 70 CONTINUE 288 ELSE 289* 290* Eliminate without interchange. 291* 292 IF( B( J, J ).EQ.ZERO ) 293 $ B( J, J ) = EPS3 294 X = CLADIV( EJ, B( J, J ) ) 295 IF( X.NE.ZERO ) THEN 296 DO 80 I = 1, J - 1 297 B( I, J-1 ) = B( I, J-1 ) - X*B( I, J ) 298 80 CONTINUE 299 END IF 300 END IF 301 90 CONTINUE 302 IF( B( 1, 1 ).EQ.ZERO ) 303 $ B( 1, 1 ) = EPS3 304* 305 TRANS = 'C' 306* 307 END IF 308* 309 NORMIN = 'N' 310 DO 110 ITS = 1, N 311* 312* Solve U*x = scale*v for a right eigenvector 313* or U**H *x = scale*v for a left eigenvector, 314* overwriting x on v. 315* 316 CALL CLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB, V, 317 $ SCALE, RWORK, IERR ) 318 NORMIN = 'Y' 319* 320* Test for sufficient growth in the norm of v. 321* 322 VNORM = SCASUM( N, V, 1 ) 323 IF( VNORM.GE.GROWTO*SCALE ) 324 $ GO TO 120 325* 326* Choose new orthogonal starting vector and try again. 327* 328 RTEMP = EPS3 / ( ROOTN+ONE ) 329 V( 1 ) = EPS3 330 DO 100 I = 2, N 331 V( I ) = RTEMP 332 100 CONTINUE 333 V( N-ITS+1 ) = V( N-ITS+1 ) - EPS3*ROOTN 334 110 CONTINUE 335* 336* Failure to find eigenvector in N iterations. 337* 338 INFO = 1 339* 340 120 CONTINUE 341* 342* Normalize eigenvector. 343* 344 I = ICAMAX( N, V, 1 ) 345 CALL CSSCAL( N, ONE / CABS1( V( I ) ), V, 1 ) 346* 347 RETURN 348* 349* End of CLAEIN 350* 351 END 352