 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B, 00002 $ LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO ) 00003 * 00004 * -- LAPACK auxiliary routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 00008 * 00009 * .. Scalar Arguments .. 00010 LOGICAL NOINIT, RIGHTV 00011 INTEGER INFO, LDB, LDH, N 00012 DOUBLE PRECISION BIGNUM, EPS3, SMLNUM, WI, WR 00013 * .. 00014 * .. Array Arguments .. 00015 DOUBLE PRECISION B( LDB, * ), H( LDH, * ), VI( * ), VR( * ), 00016 $ WORK( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * DLAEIN uses inverse iteration to find a right or left eigenvector 00023 * corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg 00024 * matrix H. 00025 * 00026 * Arguments 00027 * ========= 00028 * 00029 * RIGHTV (input) LOGICAL 00030 * = .TRUE. : compute right eigenvector; 00031 * = .FALSE.: compute left eigenvector. 00032 * 00033 * NOINIT (input) LOGICAL 00034 * = .TRUE. : no initial vector supplied in (VR,VI). 00035 * = .FALSE.: initial vector supplied in (VR,VI). 00036 * 00037 * N (input) INTEGER 00038 * The order of the matrix H. N >= 0. 00039 * 00040 * H (input) DOUBLE PRECISION array, dimension (LDH,N) 00041 * The upper Hessenberg matrix H. 00042 * 00043 * LDH (input) INTEGER 00044 * The leading dimension of the array H. LDH >= max(1,N). 00045 * 00046 * WR (input) DOUBLE PRECISION 00047 * WI (input) DOUBLE PRECISION 00048 * The real and imaginary parts of the eigenvalue of H whose 00049 * corresponding right or left eigenvector is to be computed. 00050 * 00051 * VR (input/output) DOUBLE PRECISION array, dimension (N) 00052 * VI (input/output) DOUBLE PRECISION array, dimension (N) 00053 * On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain 00054 * a real starting vector for inverse iteration using the real 00055 * eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI 00056 * must contain the real and imaginary parts of a complex 00057 * starting vector for inverse iteration using the complex 00058 * eigenvalue (WR,WI); otherwise VR and VI need not be set. 00059 * On exit, if WI = 0.0 (real eigenvalue), VR contains the 00060 * computed real eigenvector; if WI.ne.0.0 (complex eigenvalue), 00061 * VR and VI contain the real and imaginary parts of the 00062 * computed complex eigenvector. The eigenvector is normalized 00063 * so that the component of largest magnitude has magnitude 1; 00064 * here the magnitude of a complex number (x,y) is taken to be 00065 * |x| + |y|. 00066 * VI is not referenced if WI = 0.0. 00067 * 00068 * B (workspace) DOUBLE PRECISION array, dimension (LDB,N) 00069 * 00070 * LDB (input) INTEGER 00071 * The leading dimension of the array B. LDB >= N+1. 00072 * 00073 * WORK (workspace) DOUBLE PRECISION array, dimension (N) 00074 * 00075 * EPS3 (input) DOUBLE PRECISION 00076 * A small machine-dependent value which is used to perturb 00077 * close eigenvalues, and to replace zero pivots. 00078 * 00079 * SMLNUM (input) DOUBLE PRECISION 00080 * A machine-dependent value close to the underflow threshold. 00081 * 00082 * BIGNUM (input) DOUBLE PRECISION 00083 * A machine-dependent value close to the overflow threshold. 00084 * 00085 * INFO (output) INTEGER 00086 * = 0: successful exit 00087 * = 1: inverse iteration did not converge; VR is set to the 00088 * last iterate, and so is VI if WI.ne.0.0. 00089 * 00090 * ===================================================================== 00091 * 00092 * .. Parameters .. 00093 DOUBLE PRECISION ZERO, ONE, TENTH 00094 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TENTH = 1.0D-1 ) 00095 * .. 00096 * .. Local Scalars .. 00097 CHARACTER NORMIN, TRANS 00098 INTEGER I, I1, I2, I3, IERR, ITS, J 00099 DOUBLE PRECISION ABSBII, ABSBJJ, EI, EJ, GROWTO, NORM, NRMSML, 00100 $ REC, ROOTN, SCALE, TEMP, VCRIT, VMAX, VNORM, W, 00101 $ W1, X, XI, XR, Y 00102 * .. 00103 * .. External Functions .. 00104 INTEGER IDAMAX 00105 DOUBLE PRECISION DASUM, DLAPY2, DNRM2 00106 EXTERNAL IDAMAX, DASUM, DLAPY2, DNRM2 00107 * .. 00108 * .. External Subroutines .. 00109 EXTERNAL DLADIV, DLATRS, DSCAL 00110 * .. 00111 * .. Intrinsic Functions .. 00112 INTRINSIC ABS, DBLE, MAX, SQRT 00113 * .. 00114 * .. Executable Statements .. 00115 * 00116 INFO = 0 00117 * 00118 * GROWTO is the threshold used in the acceptance test for an 00119 * eigenvector. 00120 * 00121 ROOTN = SQRT( DBLE( N ) ) 00122 GROWTO = TENTH / ROOTN 00123 NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM 00124 * 00125 * Form B = H - (WR,WI)*I (except that the subdiagonal elements and 00126 * the imaginary parts of the diagonal elements are not stored). 00127 * 00128 DO 20 J = 1, N 00129 DO 10 I = 1, J - 1 00130 B( I, J ) = H( I, J ) 00131 10 CONTINUE 00132 B( J, J ) = H( J, J ) - WR 00133 20 CONTINUE 00134 * 00135 IF( WI.EQ.ZERO ) THEN 00136 * 00137 * Real eigenvalue. 00138 * 00139 IF( NOINIT ) THEN 00140 * 00141 * Set initial vector. 00142 * 00143 DO 30 I = 1, N 00144 VR( I ) = EPS3 00145 30 CONTINUE 00146 ELSE 00147 * 00148 * Scale supplied initial vector. 00149 * 00150 VNORM = DNRM2( N, VR, 1 ) 00151 CALL DSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), VR, 00152 $ 1 ) 00153 END IF 00154 * 00155 IF( RIGHTV ) THEN 00156 * 00157 * LU decomposition with partial pivoting of B, replacing zero 00158 * pivots by EPS3. 00159 * 00160 DO 60 I = 1, N - 1 00161 EI = H( I+1, I ) 00162 IF( ABS( B( I, I ) ).LT.ABS( EI ) ) THEN 00163 * 00164 * Interchange rows and eliminate. 00165 * 00166 X = B( I, I ) / EI 00167 B( I, I ) = EI 00168 DO 40 J = I + 1, N 00169 TEMP = B( I+1, J ) 00170 B( I+1, J ) = B( I, J ) - X*TEMP 00171 B( I, J ) = TEMP 00172 40 CONTINUE 00173 ELSE 00174 * 00175 * Eliminate without interchange. 00176 * 00177 IF( B( I, I ).EQ.ZERO ) 00178 $ B( I, I ) = EPS3 00179 X = EI / B( I, I ) 00180 IF( X.NE.ZERO ) THEN 00181 DO 50 J = I + 1, N 00182 B( I+1, J ) = B( I+1, J ) - X*B( I, J ) 00183 50 CONTINUE 00184 END IF 00185 END IF 00186 60 CONTINUE 00187 IF( B( N, N ).EQ.ZERO ) 00188 $ B( N, N ) = EPS3 00189 * 00190 TRANS = 'N' 00191 * 00192 ELSE 00193 * 00194 * UL decomposition with partial pivoting of B, replacing zero 00195 * pivots by EPS3. 00196 * 00197 DO 90 J = N, 2, -1 00198 EJ = H( J, J-1 ) 00199 IF( ABS( B( J, J ) ).LT.ABS( EJ ) ) THEN 00200 * 00201 * Interchange columns and eliminate. 00202 * 00203 X = B( J, J ) / EJ 00204 B( J, J ) = EJ 00205 DO 70 I = 1, J - 1 00206 TEMP = B( I, J-1 ) 00207 B( I, J-1 ) = B( I, J ) - X*TEMP 00208 B( I, J ) = TEMP 00209 70 CONTINUE 00210 ELSE 00211 * 00212 * Eliminate without interchange. 00213 * 00214 IF( B( J, J ).EQ.ZERO ) 00215 $ B( J, J ) = EPS3 00216 X = EJ / B( J, J ) 00217 IF( X.NE.ZERO ) THEN 00218 DO 80 I = 1, J - 1 00219 B( I, J-1 ) = B( I, J-1 ) - X*B( I, J ) 00220 80 CONTINUE 00221 END IF 00222 END IF 00223 90 CONTINUE 00224 IF( B( 1, 1 ).EQ.ZERO ) 00225 $ B( 1, 1 ) = EPS3 00226 * 00227 TRANS = 'T' 00228 * 00229 END IF 00230 * 00231 NORMIN = 'N' 00232 DO 110 ITS = 1, N 00233 * 00234 * Solve U*x = scale*v for a right eigenvector 00235 * or U**T*x = scale*v for a left eigenvector, 00236 * overwriting x on v. 00237 * 00238 CALL DLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB, 00239 $ VR, SCALE, WORK, IERR ) 00240 NORMIN = 'Y' 00241 * 00242 * Test for sufficient growth in the norm of v. 00243 * 00244 VNORM = DASUM( N, VR, 1 ) 00245 IF( VNORM.GE.GROWTO*SCALE ) 00246 $ GO TO 120 00247 * 00248 * Choose new orthogonal starting vector and try again. 00249 * 00250 TEMP = EPS3 / ( ROOTN+ONE ) 00251 VR( 1 ) = EPS3 00252 DO 100 I = 2, N 00253 VR( I ) = TEMP 00254 100 CONTINUE 00255 VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN 00256 110 CONTINUE 00257 * 00258 * Failure to find eigenvector in N iterations. 00259 * 00260 INFO = 1 00261 * 00262 120 CONTINUE 00263 * 00264 * Normalize eigenvector. 00265 * 00266 I = IDAMAX( N, VR, 1 ) 00267 CALL DSCAL( N, ONE / ABS( VR( I ) ), VR, 1 ) 00268 ELSE 00269 * 00270 * Complex eigenvalue. 00271 * 00272 IF( NOINIT ) THEN 00273 * 00274 * Set initial vector. 00275 * 00276 DO 130 I = 1, N 00277 VR( I ) = EPS3 00278 VI( I ) = ZERO 00279 130 CONTINUE 00280 ELSE 00281 * 00282 * Scale supplied initial vector. 00283 * 00284 NORM = DLAPY2( DNRM2( N, VR, 1 ), DNRM2( N, VI, 1 ) ) 00285 REC = ( EPS3*ROOTN ) / MAX( NORM, NRMSML ) 00286 CALL DSCAL( N, REC, VR, 1 ) 00287 CALL DSCAL( N, REC, VI, 1 ) 00288 END IF 00289 * 00290 IF( RIGHTV ) THEN 00291 * 00292 * LU decomposition with partial pivoting of B, replacing zero 00293 * pivots by EPS3. 00294 * 00295 * The imaginary part of the (i,j)-th element of U is stored in 00296 * B(j+1,i). 00297 * 00298 B( 2, 1 ) = -WI 00299 DO 140 I = 2, N 00300 B( I+1, 1 ) = ZERO 00301 140 CONTINUE 00302 * 00303 DO 170 I = 1, N - 1 00304 ABSBII = DLAPY2( B( I, I ), B( I+1, I ) ) 00305 EI = H( I+1, I ) 00306 IF( ABSBII.LT.ABS( EI ) ) THEN 00307 * 00308 * Interchange rows and eliminate. 00309 * 00310 XR = B( I, I ) / EI 00311 XI = B( I+1, I ) / EI 00312 B( I, I ) = EI 00313 B( I+1, I ) = ZERO 00314 DO 150 J = I + 1, N 00315 TEMP = B( I+1, J ) 00316 B( I+1, J ) = B( I, J ) - XR*TEMP 00317 B( J+1, I+1 ) = B( J+1, I ) - XI*TEMP 00318 B( I, J ) = TEMP 00319 B( J+1, I ) = ZERO 00320 150 CONTINUE 00321 B( I+2, I ) = -WI 00322 B( I+1, I+1 ) = B( I+1, I+1 ) - XI*WI 00323 B( I+2, I+1 ) = B( I+2, I+1 ) + XR*WI 00324 ELSE 00325 * 00326 * Eliminate without interchanging rows. 00327 * 00328 IF( ABSBII.EQ.ZERO ) THEN 00329 B( I, I ) = EPS3 00330 B( I+1, I ) = ZERO 00331 ABSBII = EPS3 00332 END IF 00333 EI = ( EI / ABSBII ) / ABSBII 00334 XR = B( I, I )*EI 00335 XI = -B( I+1, I )*EI 00336 DO 160 J = I + 1, N 00337 B( I+1, J ) = B( I+1, J ) - XR*B( I, J ) + 00338 $ XI*B( J+1, I ) 00339 B( J+1, I+1 ) = -XR*B( J+1, I ) - XI*B( I, J ) 00340 160 CONTINUE 00341 B( I+2, I+1 ) = B( I+2, I+1 ) - WI 00342 END IF 00343 * 00344 * Compute 1-norm of offdiagonal elements of i-th row. 00345 * 00346 WORK( I ) = DASUM( N-I, B( I, I+1 ), LDB ) + 00347 $ DASUM( N-I, B( I+2, I ), 1 ) 00348 170 CONTINUE 00349 IF( B( N, N ).EQ.ZERO .AND. B( N+1, N ).EQ.ZERO ) 00350 $ B( N, N ) = EPS3 00351 WORK( N ) = ZERO 00352 * 00353 I1 = N 00354 I2 = 1 00355 I3 = -1 00356 ELSE 00357 * 00358 * UL decomposition with partial pivoting of conjg(B), 00359 * replacing zero pivots by EPS3. 00360 * 00361 * The imaginary part of the (i,j)-th element of U is stored in 00362 * B(j+1,i). 00363 * 00364 B( N+1, N ) = WI 00365 DO 180 J = 1, N - 1 00366 B( N+1, J ) = ZERO 00367 180 CONTINUE 00368 * 00369 DO 210 J = N, 2, -1 00370 EJ = H( J, J-1 ) 00371 ABSBJJ = DLAPY2( B( J, J ), B( J+1, J ) ) 00372 IF( ABSBJJ.LT.ABS( EJ ) ) THEN 00373 * 00374 * Interchange columns and eliminate 00375 * 00376 XR = B( J, J ) / EJ 00377 XI = B( J+1, J ) / EJ 00378 B( J, J ) = EJ 00379 B( J+1, J ) = ZERO 00380 DO 190 I = 1, J - 1 00381 TEMP = B( I, J-1 ) 00382 B( I, J-1 ) = B( I, J ) - XR*TEMP 00383 B( J, I ) = B( J+1, I ) - XI*TEMP 00384 B( I, J ) = TEMP 00385 B( J+1, I ) = ZERO 00386 190 CONTINUE 00387 B( J+1, J-1 ) = WI 00388 B( J-1, J-1 ) = B( J-1, J-1 ) + XI*WI 00389 B( J, J-1 ) = B( J, J-1 ) - XR*WI 00390 ELSE 00391 * 00392 * Eliminate without interchange. 00393 * 00394 IF( ABSBJJ.EQ.ZERO ) THEN 00395 B( J, J ) = EPS3 00396 B( J+1, J ) = ZERO 00397 ABSBJJ = EPS3 00398 END IF 00399 EJ = ( EJ / ABSBJJ ) / ABSBJJ 00400 XR = B( J, J )*EJ 00401 XI = -B( J+1, J )*EJ 00402 DO 200 I = 1, J - 1 00403 B( I, J-1 ) = B( I, J-1 ) - XR*B( I, J ) + 00404 $ XI*B( J+1, I ) 00405 B( J, I ) = -XR*B( J+1, I ) - XI*B( I, J ) 00406 200 CONTINUE 00407 B( J, J-1 ) = B( J, J-1 ) + WI 00408 END IF 00409 * 00410 * Compute 1-norm of offdiagonal elements of j-th column. 00411 * 00412 WORK( J ) = DASUM( J-1, B( 1, J ), 1 ) + 00413 $ DASUM( J-1, B( J+1, 1 ), LDB ) 00414 210 CONTINUE 00415 IF( B( 1, 1 ).EQ.ZERO .AND. B( 2, 1 ).EQ.ZERO ) 00416 $ B( 1, 1 ) = EPS3 00417 WORK( 1 ) = ZERO 00418 * 00419 I1 = 1 00420 I2 = N 00421 I3 = 1 00422 END IF 00423 * 00424 DO 270 ITS = 1, N 00425 SCALE = ONE 00426 VMAX = ONE 00427 VCRIT = BIGNUM 00428 * 00429 * Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector, 00430 * or U**T*(xr,xi) = scale*(vr,vi) for a left eigenvector, 00431 * overwriting (xr,xi) on (vr,vi). 00432 * 00433 DO 250 I = I1, I2, I3 00434 * 00435 IF( WORK( I ).GT.VCRIT ) THEN 00436 REC = ONE / VMAX 00437 CALL DSCAL( N, REC, VR, 1 ) 00438 CALL DSCAL( N, REC, VI, 1 ) 00439 SCALE = SCALE*REC 00440 VMAX = ONE 00441 VCRIT = BIGNUM 00442 END IF 00443 * 00444 XR = VR( I ) 00445 XI = VI( I ) 00446 IF( RIGHTV ) THEN 00447 DO 220 J = I + 1, N 00448 XR = XR - B( I, J )*VR( J ) + B( J+1, I )*VI( J ) 00449 XI = XI - B( I, J )*VI( J ) - B( J+1, I )*VR( J ) 00450 220 CONTINUE 00451 ELSE 00452 DO 230 J = 1, I - 1 00453 XR = XR - B( J, I )*VR( J ) + B( I+1, J )*VI( J ) 00454 XI = XI - B( J, I )*VI( J ) - B( I+1, J )*VR( J ) 00455 230 CONTINUE 00456 END IF 00457 * 00458 W = ABS( B( I, I ) ) + ABS( B( I+1, I ) ) 00459 IF( W.GT.SMLNUM ) THEN 00460 IF( W.LT.ONE ) THEN 00461 W1 = ABS( XR ) + ABS( XI ) 00462 IF( W1.GT.W*BIGNUM ) THEN 00463 REC = ONE / W1 00464 CALL DSCAL( N, REC, VR, 1 ) 00465 CALL DSCAL( N, REC, VI, 1 ) 00466 XR = VR( I ) 00467 XI = VI( I ) 00468 SCALE = SCALE*REC 00469 VMAX = VMAX*REC 00470 END IF 00471 END IF 00472 * 00473 * Divide by diagonal element of B. 00474 * 00475 CALL DLADIV( XR, XI, B( I, I ), B( I+1, I ), VR( I ), 00476 $ VI( I ) ) 00477 VMAX = MAX( ABS( VR( I ) )+ABS( VI( I ) ), VMAX ) 00478 VCRIT = BIGNUM / VMAX 00479 ELSE 00480 DO 240 J = 1, N 00481 VR( J ) = ZERO 00482 VI( J ) = ZERO 00483 240 CONTINUE 00484 VR( I ) = ONE 00485 VI( I ) = ONE 00486 SCALE = ZERO 00487 VMAX = ONE 00488 VCRIT = BIGNUM 00489 END IF 00490 250 CONTINUE 00491 * 00492 * Test for sufficient growth in the norm of (VR,VI). 00493 * 00494 VNORM = DASUM( N, VR, 1 ) + DASUM( N, VI, 1 ) 00495 IF( VNORM.GE.GROWTO*SCALE ) 00496 $ GO TO 280 00497 * 00498 * Choose a new orthogonal starting vector and try again. 00499 * 00500 Y = EPS3 / ( ROOTN+ONE ) 00501 VR( 1 ) = EPS3 00502 VI( 1 ) = ZERO 00503 * 00504 DO 260 I = 2, N 00505 VR( I ) = Y 00506 VI( I ) = ZERO 00507 260 CONTINUE 00508 VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN 00509 270 CONTINUE 00510 * 00511 * Failure to find eigenvector in N iterations 00512 * 00513 INFO = 1 00514 * 00515 280 CONTINUE 00516 * 00517 * Normalize eigenvector. 00518 * 00519 VNORM = ZERO 00520 DO 290 I = 1, N 00521 VNORM = MAX( VNORM, ABS( VR( I ) )+ABS( VI( I ) ) ) 00522 290 CONTINUE 00523 CALL DSCAL( N, ONE / VNORM, VR, 1 ) 00524 CALL DSCAL( N, ONE / VNORM, VI, 1 ) 00525 * 00526 END IF 00527 * 00528 RETURN 00529 * 00530 * End of DLAEIN 00531 * 00532 END
1.7.3 
