|
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 |
SUBROUTINE ZHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL,
$ LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL, $ IFAILR, INFO ) * * -- LAPACK routine (version 3.2) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2006 * * .. Scalar Arguments .. CHARACTER EIGSRC, INITV, SIDE INTEGER INFO, LDH, LDVL, LDVR, M, MM, N * .. * .. Array Arguments .. LOGICAL SELECT( * ) INTEGER IFAILL( * ), IFAILR( * ) DOUBLE PRECISION RWORK( * ) COMPLEX*16 H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ), $ W( * ), WORK( * ) * .. * * Purpose * ======= * * ZHSEIN uses inverse iteration to find specified right and/or left * eigenvectors of a complex upper Hessenberg matrix H. * * The right eigenvector x and the left eigenvector y of the matrix H * corresponding to an eigenvalue w are defined by: * * H * x = w * x, y**h * H = w * y**h * * where y**h denotes the conjugate transpose of the vector y. * * Arguments * ========= * * SIDE (input) CHARACTER*1 * = 'R': compute right eigenvectors only; * = 'L': compute left eigenvectors only; * = 'B': compute both right and left eigenvectors. * * EIGSRC (input) CHARACTER*1 * Specifies the source of eigenvalues supplied in W: * = 'Q': the eigenvalues were found using ZHSEQR; thus, if * H has zero subdiagonal elements, and so is * block-triangular, then the j-th eigenvalue can be * assumed to be an eigenvalue of the block containing * the j-th row/column. This property allows ZHSEIN to * perform inverse iteration on just one diagonal block. * = 'N': no assumptions are made on the correspondence * between eigenvalues and diagonal blocks. In this * case, ZHSEIN must always perform inverse iteration * using the whole matrix H. * * INITV (input) CHARACTER*1 * = 'N': no initial vectors are supplied; * = 'U': user-supplied initial vectors are stored in the arrays * VL and/or VR. * * SELECT (input) LOGICAL array, dimension (N) * Specifies the eigenvectors to be computed. To select the * eigenvector corresponding to the eigenvalue W(j), * SELECT(j) must be set to .TRUE.. * * N (input) INTEGER * The order of the matrix H. N >= 0. * * H (input) COMPLEX*16 array, dimension (LDH,N) * The upper Hessenberg matrix H. * * LDH (input) INTEGER * The leading dimension of the array H. LDH >= max(1,N). * * W (input/output) COMPLEX*16 array, dimension (N) * On entry, the eigenvalues of H. * On exit, the real parts of W may have been altered since * close eigenvalues are perturbed slightly in searching for * independent eigenvectors. * * VL (input/output) COMPLEX*16 array, dimension (LDVL,MM) * On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must * contain starting vectors for the inverse iteration for the * left eigenvectors; the starting vector for each eigenvector * must be in the same column in which the eigenvector will be * stored. * On exit, if SIDE = 'L' or 'B', the left eigenvectors * specified by SELECT will be stored consecutively in the * columns of VL, in the same order as their eigenvalues. * If SIDE = 'R', VL is not referenced. * * LDVL (input) INTEGER * The leading dimension of the array VL. * LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise. * * VR (input/output) COMPLEX*16 array, dimension (LDVR,MM) * On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must * contain starting vectors for the inverse iteration for the * right eigenvectors; the starting vector for each eigenvector * must be in the same column in which the eigenvector will be * stored. * On exit, if SIDE = 'R' or 'B', the right eigenvectors * specified by SELECT will be stored consecutively in the * columns of VR, in the same order as their eigenvalues. * If SIDE = 'L', VR is not referenced. * * LDVR (input) INTEGER * The leading dimension of the array VR. * LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise. * * MM (input) INTEGER * The number of columns in the arrays VL and/or VR. MM >= M. * * M (output) INTEGER * The number of columns in the arrays VL and/or VR required to * store the eigenvectors (= the number of .TRUE. elements in * SELECT). * * WORK (workspace) COMPLEX*16 array, dimension (N*N) * * RWORK (workspace) DOUBLE PRECISION array, dimension (N) * * IFAILL (output) INTEGER array, dimension (MM) * If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left * eigenvector in the i-th column of VL (corresponding to the * eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the * eigenvector converged satisfactorily. * If SIDE = 'R', IFAILL is not referenced. * * IFAILR (output) INTEGER array, dimension (MM) * If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right * eigenvector in the i-th column of VR (corresponding to the * eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the * eigenvector converged satisfactorily. * If SIDE = 'L', IFAILR is not referenced. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: if INFO = i, i is the number of eigenvectors which * failed to converge; see IFAILL and IFAILR for further * details. * * Further Details * =============== * * Each eigenvector is normalized so that the element of largest * magnitude has magnitude 1; here the magnitude of a complex number * (x,y) is taken to be |x|+|y|. * * ===================================================================== * * .. Parameters .. COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) DOUBLE PRECISION RZERO PARAMETER ( RZERO = 0.0D+0 ) * .. * .. Local Scalars .. LOGICAL BOTHV, FROMQR, LEFTV, NOINIT, RIGHTV INTEGER I, IINFO, K, KL, KLN, KR, KS, LDWORK DOUBLE PRECISION EPS3, HNORM, SMLNUM, ULP, UNFL COMPLEX*16 CDUM, WK * .. * .. External Functions .. LOGICAL LSAME DOUBLE PRECISION DLAMCH, ZLANHS EXTERNAL LSAME, DLAMCH, ZLANHS * .. * .. External Subroutines .. EXTERNAL XERBLA, ZLAEIN * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, DIMAG, MAX * .. * .. Statement Functions .. DOUBLE PRECISION CABS1 * .. * .. Statement Function definitions .. CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) * .. * .. Executable Statements .. * * Decode and test the input parameters. * BOTHV = LSAME( SIDE, 'B' ) RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV * FROMQR = LSAME( EIGSRC, 'Q' ) * NOINIT = LSAME( INITV, 'N' ) * * Set M to the number of columns required to store the selected * eigenvectors. * M = 0 DO 10 K = 1, N IF( SELECT( K ) ) $ M = M + 1 10 CONTINUE * INFO = 0 IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN INFO = -1 ELSE IF( .NOT.FROMQR .AND. .NOT.LSAME( EIGSRC, 'N' ) ) THEN INFO = -2 ELSE IF( .NOT.NOINIT .AND. .NOT.LSAME( INITV, 'U' ) ) THEN INFO = -3 ELSE IF( N.LT.0 ) THEN INFO = -5 ELSE IF( LDH.LT.MAX( 1, N ) ) THEN INFO = -7 ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN INFO = -10 ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN INFO = -12 ELSE IF( MM.LT.M ) THEN INFO = -13 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'ZHSEIN', -INFO ) RETURN END IF * * Quick return if possible. * IF( N.EQ.0 ) $ RETURN * * Set machine-dependent constants. * UNFL = DLAMCH( 'Safe minimum' ) ULP = DLAMCH( 'Precision' ) SMLNUM = UNFL*( N / ULP ) * LDWORK = N * KL = 1 KLN = 0 IF( FROMQR ) THEN KR = 0 ELSE KR = N END IF KS = 1 * DO 100 K = 1, N IF( SELECT( K ) ) THEN * * Compute eigenvector(s) corresponding to W(K). * IF( FROMQR ) THEN * * If affiliation of eigenvalues is known, check whether * the matrix splits. * * Determine KL and KR such that 1 <= KL <= K <= KR <= N * and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or * KR = N). * * Then inverse iteration can be performed with the * submatrix H(KL:N,KL:N) for a left eigenvector, and with * the submatrix H(1:KR,1:KR) for a right eigenvector. * DO 20 I = K, KL + 1, -1 IF( H( I, I-1 ).EQ.ZERO ) $ GO TO 30 20 CONTINUE 30 CONTINUE KL = I IF( K.GT.KR ) THEN DO 40 I = K, N - 1 IF( H( I+1, I ).EQ.ZERO ) $ GO TO 50 40 CONTINUE 50 CONTINUE KR = I END IF END IF * IF( KL.NE.KLN ) THEN KLN = KL * * Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it * has not ben computed before. * HNORM = ZLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, RWORK ) IF( HNORM.GT.RZERO ) THEN EPS3 = HNORM*ULP ELSE EPS3 = SMLNUM END IF END IF * * Perturb eigenvalue if it is close to any previous * selected eigenvalues affiliated to the submatrix * H(KL:KR,KL:KR). Close roots are modified by EPS3. * WK = W( K ) 60 CONTINUE DO 70 I = K - 1, KL, -1 IF( SELECT( I ) .AND. CABS1( W( I )-WK ).LT.EPS3 ) THEN WK = WK + EPS3 GO TO 60 END IF 70 CONTINUE W( K ) = WK * IF( LEFTV ) THEN * * Compute left eigenvector. * CALL ZLAEIN( .FALSE., NOINIT, N-KL+1, H( KL, KL ), LDH, $ WK, VL( KL, KS ), WORK, LDWORK, RWORK, EPS3, $ SMLNUM, IINFO ) IF( IINFO.GT.0 ) THEN INFO = INFO + 1 IFAILL( KS ) = K ELSE IFAILL( KS ) = 0 END IF DO 80 I = 1, KL - 1 VL( I, KS ) = ZERO 80 CONTINUE END IF IF( RIGHTV ) THEN * * Compute right eigenvector. * CALL ZLAEIN( .TRUE., NOINIT, KR, H, LDH, WK, VR( 1, KS ), $ WORK, LDWORK, RWORK, EPS3, SMLNUM, IINFO ) IF( IINFO.GT.0 ) THEN INFO = INFO + 1 IFAILR( KS ) = K ELSE IFAILR( KS ) = 0 END IF DO 90 I = KR + 1, N VR( I, KS ) = ZERO 90 CONTINUE END IF KS = KS + 1 END IF 100 CONTINUE * RETURN * * End of ZHSEIN * END |