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