00001 SUBROUTINE SHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, WR, WI, 00002 $ VL, LDVL, VR, LDVR, MM, M, WORK, IFAILL, 00003 $ IFAILR, INFO ) 00004 * 00005 * -- LAPACK routine (version 3.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * November 2006 00009 * 00010 * .. Scalar Arguments .. 00011 CHARACTER EIGSRC, INITV, SIDE 00012 INTEGER INFO, LDH, LDVL, LDVR, M, MM, N 00013 * .. 00014 * .. Array Arguments .. 00015 LOGICAL SELECT( * ) 00016 INTEGER IFAILL( * ), IFAILR( * ) 00017 REAL H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ), 00018 $ WI( * ), WORK( * ), WR( * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * SHSEIN uses inverse iteration to find specified right and/or left 00025 * eigenvectors of a real upper Hessenberg matrix H. 00026 * 00027 * The right eigenvector x and the left eigenvector y of the matrix H 00028 * corresponding to an eigenvalue w are defined by: 00029 * 00030 * H * x = w * x, y**h * H = w * y**h 00031 * 00032 * where y**h denotes the conjugate transpose of the vector y. 00033 * 00034 * Arguments 00035 * ========= 00036 * 00037 * SIDE (input) CHARACTER*1 00038 * = 'R': compute right eigenvectors only; 00039 * = 'L': compute left eigenvectors only; 00040 * = 'B': compute both right and left eigenvectors. 00041 * 00042 * EIGSRC (input) CHARACTER*1 00043 * Specifies the source of eigenvalues supplied in (WR,WI): 00044 * = 'Q': the eigenvalues were found using SHSEQR; thus, if 00045 * H has zero subdiagonal elements, and so is 00046 * block-triangular, then the j-th eigenvalue can be 00047 * assumed to be an eigenvalue of the block containing 00048 * the j-th row/column. This property allows SHSEIN to 00049 * perform inverse iteration on just one diagonal block. 00050 * = 'N': no assumptions are made on the correspondence 00051 * between eigenvalues and diagonal blocks. In this 00052 * case, SHSEIN must always perform inverse iteration 00053 * using the whole matrix H. 00054 * 00055 * INITV (input) CHARACTER*1 00056 * = 'N': no initial vectors are supplied; 00057 * = 'U': user-supplied initial vectors are stored in the arrays 00058 * VL and/or VR. 00059 * 00060 * SELECT (input/output) LOGICAL array, dimension (N) 00061 * Specifies the eigenvectors to be computed. To select the 00062 * real eigenvector corresponding to a real eigenvalue WR(j), 00063 * SELECT(j) must be set to .TRUE.. To select the complex 00064 * eigenvector corresponding to a complex eigenvalue 00065 * (WR(j),WI(j)), with complex conjugate (WR(j+1),WI(j+1)), 00066 * either SELECT(j) or SELECT(j+1) or both must be set to 00067 * .TRUE.; then on exit SELECT(j) is .TRUE. and SELECT(j+1) is 00068 * .FALSE.. 00069 * 00070 * N (input) INTEGER 00071 * The order of the matrix H. N >= 0. 00072 * 00073 * H (input) REAL array, dimension (LDH,N) 00074 * The upper Hessenberg matrix H. 00075 * 00076 * LDH (input) INTEGER 00077 * The leading dimension of the array H. LDH >= max(1,N). 00078 * 00079 * WR (input/output) REAL array, dimension (N) 00080 * WI (input) REAL array, dimension (N) 00081 * On entry, the real and imaginary parts of the eigenvalues of 00082 * H; a complex conjugate pair of eigenvalues must be stored in 00083 * consecutive elements of WR and WI. 00084 * On exit, WR may have been altered since close eigenvalues 00085 * are perturbed slightly in searching for independent 00086 * eigenvectors. 00087 * 00088 * VL (input/output) REAL array, dimension (LDVL,MM) 00089 * On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must 00090 * contain starting vectors for the inverse iteration for the 00091 * left eigenvectors; the starting vector for each eigenvector 00092 * must be in the same column(s) in which the eigenvector will 00093 * be stored. 00094 * On exit, if SIDE = 'L' or 'B', the left eigenvectors 00095 * specified by SELECT will be stored consecutively in the 00096 * columns of VL, in the same order as their eigenvalues. A 00097 * complex eigenvector corresponding to a complex eigenvalue is 00098 * stored in two consecutive columns, the first holding the real 00099 * part and the second the imaginary part. 00100 * If SIDE = 'R', VL is not referenced. 00101 * 00102 * LDVL (input) INTEGER 00103 * The leading dimension of the array VL. 00104 * LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise. 00105 * 00106 * VR (input/output) REAL array, dimension (LDVR,MM) 00107 * On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must 00108 * contain starting vectors for the inverse iteration for the 00109 * right eigenvectors; the starting vector for each eigenvector 00110 * must be in the same column(s) in which the eigenvector will 00111 * be stored. 00112 * On exit, if SIDE = 'R' or 'B', the right eigenvectors 00113 * specified by SELECT will be stored consecutively in the 00114 * columns of VR, in the same order as their eigenvalues. A 00115 * complex eigenvector corresponding to a complex eigenvalue is 00116 * stored in two consecutive columns, the first holding the real 00117 * part and the second the imaginary part. 00118 * If SIDE = 'L', VR is not referenced. 00119 * 00120 * LDVR (input) INTEGER 00121 * The leading dimension of the array VR. 00122 * LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise. 00123 * 00124 * MM (input) INTEGER 00125 * The number of columns in the arrays VL and/or VR. MM >= M. 00126 * 00127 * M (output) INTEGER 00128 * The number of columns in the arrays VL and/or VR required to 00129 * store the eigenvectors; each selected real eigenvector 00130 * occupies one column and each selected complex eigenvector 00131 * occupies two columns. 00132 * 00133 * WORK (workspace) REAL array, dimension ((N+2)*N) 00134 * 00135 * IFAILL (output) INTEGER array, dimension (MM) 00136 * If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left 00137 * eigenvector in the i-th column of VL (corresponding to the 00138 * eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the 00139 * eigenvector converged satisfactorily. If the i-th and (i+1)th 00140 * columns of VL hold a complex eigenvector, then IFAILL(i) and 00141 * IFAILL(i+1) are set to the same value. 00142 * If SIDE = 'R', IFAILL is not referenced. 00143 * 00144 * IFAILR (output) INTEGER array, dimension (MM) 00145 * If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right 00146 * eigenvector in the i-th column of VR (corresponding to the 00147 * eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the 00148 * eigenvector converged satisfactorily. If the i-th and (i+1)th 00149 * columns of VR hold a complex eigenvector, then IFAILR(i) and 00150 * IFAILR(i+1) are set to the same value. 00151 * If SIDE = 'L', IFAILR is not referenced. 00152 * 00153 * INFO (output) INTEGER 00154 * = 0: successful exit 00155 * < 0: if INFO = -i, the i-th argument had an illegal value 00156 * > 0: if INFO = i, i is the number of eigenvectors which 00157 * failed to converge; see IFAILL and IFAILR for further 00158 * details. 00159 * 00160 * Further Details 00161 * =============== 00162 * 00163 * Each eigenvector is normalized so that the element of largest 00164 * magnitude has magnitude 1; here the magnitude of a complex number 00165 * (x,y) is taken to be |x|+|y|. 00166 * 00167 * ===================================================================== 00168 * 00169 * .. Parameters .. 00170 REAL ZERO, ONE 00171 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00172 * .. 00173 * .. Local Scalars .. 00174 LOGICAL BOTHV, FROMQR, LEFTV, NOINIT, PAIR, RIGHTV 00175 INTEGER I, IINFO, K, KL, KLN, KR, KSI, KSR, LDWORK 00176 REAL BIGNUM, EPS3, HNORM, SMLNUM, ULP, UNFL, WKI, 00177 $ WKR 00178 * .. 00179 * .. External Functions .. 00180 LOGICAL LSAME 00181 REAL SLAMCH, SLANHS 00182 EXTERNAL LSAME, SLAMCH, SLANHS 00183 * .. 00184 * .. External Subroutines .. 00185 EXTERNAL SLAEIN, XERBLA 00186 * .. 00187 * .. Intrinsic Functions .. 00188 INTRINSIC ABS, MAX 00189 * .. 00190 * .. Executable Statements .. 00191 * 00192 * Decode and test the input parameters. 00193 * 00194 BOTHV = LSAME( SIDE, 'B' ) 00195 RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV 00196 LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV 00197 * 00198 FROMQR = LSAME( EIGSRC, 'Q' ) 00199 * 00200 NOINIT = LSAME( INITV, 'N' ) 00201 * 00202 * Set M to the number of columns required to store the selected 00203 * eigenvectors, and standardize the array SELECT. 00204 * 00205 M = 0 00206 PAIR = .FALSE. 00207 DO 10 K = 1, N 00208 IF( PAIR ) THEN 00209 PAIR = .FALSE. 00210 SELECT( K ) = .FALSE. 00211 ELSE 00212 IF( WI( K ).EQ.ZERO ) THEN 00213 IF( SELECT( K ) ) 00214 $ M = M + 1 00215 ELSE 00216 PAIR = .TRUE. 00217 IF( SELECT( K ) .OR. SELECT( K+1 ) ) THEN 00218 SELECT( K ) = .TRUE. 00219 M = M + 2 00220 END IF 00221 END IF 00222 END IF 00223 10 CONTINUE 00224 * 00225 INFO = 0 00226 IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN 00227 INFO = -1 00228 ELSE IF( .NOT.FROMQR .AND. .NOT.LSAME( EIGSRC, 'N' ) ) THEN 00229 INFO = -2 00230 ELSE IF( .NOT.NOINIT .AND. .NOT.LSAME( INITV, 'U' ) ) THEN 00231 INFO = -3 00232 ELSE IF( N.LT.0 ) THEN 00233 INFO = -5 00234 ELSE IF( LDH.LT.MAX( 1, N ) ) THEN 00235 INFO = -7 00236 ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN 00237 INFO = -11 00238 ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN 00239 INFO = -13 00240 ELSE IF( MM.LT.M ) THEN 00241 INFO = -14 00242 END IF 00243 IF( INFO.NE.0 ) THEN 00244 CALL XERBLA( 'SHSEIN', -INFO ) 00245 RETURN 00246 END IF 00247 * 00248 * Quick return if possible. 00249 * 00250 IF( N.EQ.0 ) 00251 $ RETURN 00252 * 00253 * Set machine-dependent constants. 00254 * 00255 UNFL = SLAMCH( 'Safe minimum' ) 00256 ULP = SLAMCH( 'Precision' ) 00257 SMLNUM = UNFL*( N / ULP ) 00258 BIGNUM = ( ONE-ULP ) / SMLNUM 00259 * 00260 LDWORK = N + 1 00261 * 00262 KL = 1 00263 KLN = 0 00264 IF( FROMQR ) THEN 00265 KR = 0 00266 ELSE 00267 KR = N 00268 END IF 00269 KSR = 1 00270 * 00271 DO 120 K = 1, N 00272 IF( SELECT( K ) ) THEN 00273 * 00274 * Compute eigenvector(s) corresponding to W(K). 00275 * 00276 IF( FROMQR ) THEN 00277 * 00278 * If affiliation of eigenvalues is known, check whether 00279 * the matrix splits. 00280 * 00281 * Determine KL and KR such that 1 <= KL <= K <= KR <= N 00282 * and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or 00283 * KR = N). 00284 * 00285 * Then inverse iteration can be performed with the 00286 * submatrix H(KL:N,KL:N) for a left eigenvector, and with 00287 * the submatrix H(1:KR,1:KR) for a right eigenvector. 00288 * 00289 DO 20 I = K, KL + 1, -1 00290 IF( H( I, I-1 ).EQ.ZERO ) 00291 $ GO TO 30 00292 20 CONTINUE 00293 30 CONTINUE 00294 KL = I 00295 IF( K.GT.KR ) THEN 00296 DO 40 I = K, N - 1 00297 IF( H( I+1, I ).EQ.ZERO ) 00298 $ GO TO 50 00299 40 CONTINUE 00300 50 CONTINUE 00301 KR = I 00302 END IF 00303 END IF 00304 * 00305 IF( KL.NE.KLN ) THEN 00306 KLN = KL 00307 * 00308 * Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it 00309 * has not ben computed before. 00310 * 00311 HNORM = SLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, WORK ) 00312 IF( HNORM.GT.ZERO ) THEN 00313 EPS3 = HNORM*ULP 00314 ELSE 00315 EPS3 = SMLNUM 00316 END IF 00317 END IF 00318 * 00319 * Perturb eigenvalue if it is close to any previous 00320 * selected eigenvalues affiliated to the submatrix 00321 * H(KL:KR,KL:KR). Close roots are modified by EPS3. 00322 * 00323 WKR = WR( K ) 00324 WKI = WI( K ) 00325 60 CONTINUE 00326 DO 70 I = K - 1, KL, -1 00327 IF( SELECT( I ) .AND. ABS( WR( I )-WKR )+ 00328 $ ABS( WI( I )-WKI ).LT.EPS3 ) THEN 00329 WKR = WKR + EPS3 00330 GO TO 60 00331 END IF 00332 70 CONTINUE 00333 WR( K ) = WKR 00334 * 00335 PAIR = WKI.NE.ZERO 00336 IF( PAIR ) THEN 00337 KSI = KSR + 1 00338 ELSE 00339 KSI = KSR 00340 END IF 00341 IF( LEFTV ) THEN 00342 * 00343 * Compute left eigenvector. 00344 * 00345 CALL SLAEIN( .FALSE., NOINIT, N-KL+1, H( KL, KL ), LDH, 00346 $ WKR, WKI, VL( KL, KSR ), VL( KL, KSI ), 00347 $ WORK, LDWORK, WORK( N*N+N+1 ), EPS3, SMLNUM, 00348 $ BIGNUM, IINFO ) 00349 IF( IINFO.GT.0 ) THEN 00350 IF( PAIR ) THEN 00351 INFO = INFO + 2 00352 ELSE 00353 INFO = INFO + 1 00354 END IF 00355 IFAILL( KSR ) = K 00356 IFAILL( KSI ) = K 00357 ELSE 00358 IFAILL( KSR ) = 0 00359 IFAILL( KSI ) = 0 00360 END IF 00361 DO 80 I = 1, KL - 1 00362 VL( I, KSR ) = ZERO 00363 80 CONTINUE 00364 IF( PAIR ) THEN 00365 DO 90 I = 1, KL - 1 00366 VL( I, KSI ) = ZERO 00367 90 CONTINUE 00368 END IF 00369 END IF 00370 IF( RIGHTV ) THEN 00371 * 00372 * Compute right eigenvector. 00373 * 00374 CALL SLAEIN( .TRUE., NOINIT, KR, H, LDH, WKR, WKI, 00375 $ VR( 1, KSR ), VR( 1, KSI ), WORK, LDWORK, 00376 $ WORK( N*N+N+1 ), EPS3, SMLNUM, BIGNUM, 00377 $ IINFO ) 00378 IF( IINFO.GT.0 ) THEN 00379 IF( PAIR ) THEN 00380 INFO = INFO + 2 00381 ELSE 00382 INFO = INFO + 1 00383 END IF 00384 IFAILR( KSR ) = K 00385 IFAILR( KSI ) = K 00386 ELSE 00387 IFAILR( KSR ) = 0 00388 IFAILR( KSI ) = 0 00389 END IF 00390 DO 100 I = KR + 1, N 00391 VR( I, KSR ) = ZERO 00392 100 CONTINUE 00393 IF( PAIR ) THEN 00394 DO 110 I = KR + 1, N 00395 VR( I, KSI ) = ZERO 00396 110 CONTINUE 00397 END IF 00398 END IF 00399 * 00400 IF( PAIR ) THEN 00401 KSR = KSR + 2 00402 ELSE 00403 KSR = KSR + 1 00404 END IF 00405 END IF 00406 120 CONTINUE 00407 * 00408 RETURN 00409 * 00410 * End of SHSEIN 00411 * 00412 END
1.7.2