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