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