![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZHSEIN 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZHSEIN + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhsein.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhsein.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhsein.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZHSEIN( 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 * DOUBLE PRECISION RWORK( * ) 00033 * COMPLEX*16 H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ), 00034 * $ W( * ), WORK( * ) 00035 * .. 00036 * 00037 * 00038 *> \par Purpose: 00039 * ============= 00040 *> 00041 *> \verbatim 00042 *> 00043 *> ZHSEIN 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 ZHSEQR; 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 ZHSEIN 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, ZHSEIN 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*16 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*16 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*16 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*16 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*16 array, dimension (N*N) 00183 *> \endverbatim 00184 *> 00185 *> \param[out] RWORK 00186 *> \verbatim 00187 *> RWORK is DOUBLE PRECISION 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 complex16OTHERcomputational 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 ZHSEIN( 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 DOUBLE PRECISION RWORK( * ) 00260 COMPLEX*16 H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ), 00261 $ W( * ), WORK( * ) 00262 * .. 00263 * 00264 * ===================================================================== 00265 * 00266 * .. Parameters .. 00267 COMPLEX*16 ZERO 00268 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) 00269 DOUBLE PRECISION RZERO 00270 PARAMETER ( RZERO = 0.0D+0 ) 00271 * .. 00272 * .. Local Scalars .. 00273 LOGICAL BOTHV, FROMQR, LEFTV, NOINIT, RIGHTV 00274 INTEGER I, IINFO, K, KL, KLN, KR, KS, LDWORK 00275 DOUBLE PRECISION EPS3, HNORM, SMLNUM, ULP, UNFL 00276 COMPLEX*16 CDUM, WK 00277 * .. 00278 * .. External Functions .. 00279 LOGICAL LSAME 00280 DOUBLE PRECISION DLAMCH, ZLANHS 00281 EXTERNAL LSAME, DLAMCH, ZLANHS 00282 * .. 00283 * .. External Subroutines .. 00284 EXTERNAL XERBLA, ZLAEIN 00285 * .. 00286 * .. Intrinsic Functions .. 00287 INTRINSIC ABS, DBLE, DIMAG, MAX 00288 * .. 00289 * .. Statement Functions .. 00290 DOUBLE PRECISION CABS1 00291 * .. 00292 * .. Statement Function definitions .. 00293 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( 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( 'ZHSEIN', -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 = DLAMCH( 'Safe minimum' ) 00347 ULP = DLAMCH( '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 = ZLANHS( '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 ZLAEIN( .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 ZLAEIN( .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 ZHSEIN 00463 * 00464 END