 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CHSEIN( 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 REAL RWORK( * ) 00018 COMPLEX H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ), 00019 $ W( * ), WORK( * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * CHSEIN 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 CHSEQR; 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 CHSEIN 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, CHSEIN 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 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 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 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 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 array, dimension (N*N) 00120 * 00121 * RWORK (workspace) REAL 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 ZERO 00155 PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ) ) 00156 REAL RZERO 00157 PARAMETER ( RZERO = 0.0E+0 ) 00158 * .. 00159 * .. Local Scalars .. 00160 LOGICAL BOTHV, FROMQR, LEFTV, NOINIT, RIGHTV 00161 INTEGER I, IINFO, K, KL, KLN, KR, KS, LDWORK 00162 REAL EPS3, HNORM, SMLNUM, ULP, UNFL 00163 COMPLEX CDUM, WK 00164 * .. 00165 * .. External Functions .. 00166 LOGICAL LSAME 00167 REAL CLANHS, SLAMCH 00168 EXTERNAL LSAME, CLANHS, SLAMCH 00169 * .. 00170 * .. External Subroutines .. 00171 EXTERNAL CLAEIN, XERBLA 00172 * .. 00173 * .. Intrinsic Functions .. 00174 INTRINSIC ABS, AIMAG, MAX, REAL 00175 * .. 00176 * .. Statement Functions .. 00177 REAL CABS1 00178 * .. 00179 * .. Statement Function definitions .. 00180 CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( 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( 'CHSEIN', -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 = SLAMCH( 'Safe minimum' ) 00234 ULP = SLAMCH( '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 = CLANHS( '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 CLAEIN( .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 CLAEIN( .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 CHSEIN 00350 * 00351 END
1.7.3 
