001: SUBROUTINE ZHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL, 002: $ LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL, 003: $ IFAILR, INFO ) 004: * 005: * -- LAPACK routine (version 3.2) -- 006: * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 007: * November 2006 008: * 009: * .. Scalar Arguments .. 010: CHARACTER EIGSRC, INITV, SIDE 011: INTEGER INFO, LDH, LDVL, LDVR, M, MM, N 012: * .. 013: * .. Array Arguments .. 014: LOGICAL SELECT( * ) 015: INTEGER IFAILL( * ), IFAILR( * ) 016: DOUBLE PRECISION RWORK( * ) 017: COMPLEX*16 H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ), 018: $ W( * ), WORK( * ) 019: * .. 020: * 021: * Purpose 022: * ======= 023: * 024: * ZHSEIN uses inverse iteration to find specified right and/or left 025: * eigenvectors of a complex upper Hessenberg matrix H. 026: * 027: * The right eigenvector x and the left eigenvector y of the matrix H 028: * corresponding to an eigenvalue w are defined by: 029: * 030: * H * x = w * x, y**h * H = w * y**h 031: * 032: * where y**h denotes the conjugate transpose of the vector y. 033: * 034: * Arguments 035: * ========= 036: * 037: * SIDE (input) CHARACTER*1 038: * = 'R': compute right eigenvectors only; 039: * = 'L': compute left eigenvectors only; 040: * = 'B': compute both right and left eigenvectors. 041: * 042: * EIGSRC (input) CHARACTER*1 043: * Specifies the source of eigenvalues supplied in W: 044: * = 'Q': the eigenvalues were found using ZHSEQR; thus, if 045: * H has zero subdiagonal elements, and so is 046: * block-triangular, then the j-th eigenvalue can be 047: * assumed to be an eigenvalue of the block containing 048: * the j-th row/column. This property allows ZHSEIN to 049: * perform inverse iteration on just one diagonal block. 050: * = 'N': no assumptions are made on the correspondence 051: * between eigenvalues and diagonal blocks. In this 052: * case, ZHSEIN must always perform inverse iteration 053: * using the whole matrix H. 054: * 055: * INITV (input) CHARACTER*1 056: * = 'N': no initial vectors are supplied; 057: * = 'U': user-supplied initial vectors are stored in the arrays 058: * VL and/or VR. 059: * 060: * SELECT (input) LOGICAL array, dimension (N) 061: * Specifies the eigenvectors to be computed. To select the 062: * eigenvector corresponding to the eigenvalue W(j), 063: * SELECT(j) must be set to .TRUE.. 064: * 065: * N (input) INTEGER 066: * The order of the matrix H. N >= 0. 067: * 068: * H (input) COMPLEX*16 array, dimension (LDH,N) 069: * The upper Hessenberg matrix H. 070: * 071: * LDH (input) INTEGER 072: * The leading dimension of the array H. LDH >= max(1,N). 073: * 074: * W (input/output) COMPLEX*16 array, dimension (N) 075: * On entry, the eigenvalues of H. 076: * On exit, the real parts of W may have been altered since 077: * close eigenvalues are perturbed slightly in searching for 078: * independent eigenvectors. 079: * 080: * VL (input/output) COMPLEX*16 array, dimension (LDVL,MM) 081: * On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must 082: * contain starting vectors for the inverse iteration for the 083: * left eigenvectors; the starting vector for each eigenvector 084: * must be in the same column in which the eigenvector will be 085: * stored. 086: * On exit, if SIDE = 'L' or 'B', the left eigenvectors 087: * specified by SELECT will be stored consecutively in the 088: * columns of VL, in the same order as their eigenvalues. 089: * If SIDE = 'R', VL is not referenced. 090: * 091: * LDVL (input) INTEGER 092: * The leading dimension of the array VL. 093: * LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise. 094: * 095: * VR (input/output) COMPLEX*16 array, dimension (LDVR,MM) 096: * On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must 097: * contain starting vectors for the inverse iteration for the 098: * right eigenvectors; the starting vector for each eigenvector 099: * must be in the same column in which the eigenvector will be 100: * stored. 101: * On exit, if SIDE = 'R' or 'B', the right eigenvectors 102: * specified by SELECT will be stored consecutively in the 103: * columns of VR, in the same order as their eigenvalues. 104: * If SIDE = 'L', VR is not referenced. 105: * 106: * LDVR (input) INTEGER 107: * The leading dimension of the array VR. 108: * LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise. 109: * 110: * MM (input) INTEGER 111: * The number of columns in the arrays VL and/or VR. MM >= M. 112: * 113: * M (output) INTEGER 114: * The number of columns in the arrays VL and/or VR required to 115: * store the eigenvectors (= the number of .TRUE. elements in 116: * SELECT). 117: * 118: * WORK (workspace) COMPLEX*16 array, dimension (N*N) 119: * 120: * RWORK (workspace) DOUBLE PRECISION array, dimension (N) 121: * 122: * IFAILL (output) INTEGER array, dimension (MM) 123: * If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left 124: * eigenvector in the i-th column of VL (corresponding to the 125: * eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the 126: * eigenvector converged satisfactorily. 127: * If SIDE = 'R', IFAILL is not referenced. 128: * 129: * IFAILR (output) INTEGER array, dimension (MM) 130: * If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right 131: * eigenvector in the i-th column of VR (corresponding to the 132: * eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the 133: * eigenvector converged satisfactorily. 134: * If SIDE = 'L', IFAILR is not referenced. 135: * 136: * INFO (output) INTEGER 137: * = 0: successful exit 138: * < 0: if INFO = -i, the i-th argument had an illegal value 139: * > 0: if INFO = i, i is the number of eigenvectors which 140: * failed to converge; see IFAILL and IFAILR for further 141: * details. 142: * 143: * Further Details 144: * =============== 145: * 146: * Each eigenvector is normalized so that the element of largest 147: * magnitude has magnitude 1; here the magnitude of a complex number 148: * (x,y) is taken to be |x|+|y|. 149: * 150: * ===================================================================== 151: * 152: * .. Parameters .. 153: COMPLEX*16 ZERO 154: PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) 155: DOUBLE PRECISION RZERO 156: PARAMETER ( RZERO = 0.0D+0 ) 157: * .. 158: * .. Local Scalars .. 159: LOGICAL BOTHV, FROMQR, LEFTV, NOINIT, RIGHTV 160: INTEGER I, IINFO, K, KL, KLN, KR, KS, LDWORK 161: DOUBLE PRECISION EPS3, HNORM, SMLNUM, ULP, UNFL 162: COMPLEX*16 CDUM, WK 163: * .. 164: * .. External Functions .. 165: LOGICAL LSAME 166: DOUBLE PRECISION DLAMCH, ZLANHS 167: EXTERNAL LSAME, DLAMCH, ZLANHS 168: * .. 169: * .. External Subroutines .. 170: EXTERNAL XERBLA, ZLAEIN 171: * .. 172: * .. Intrinsic Functions .. 173: INTRINSIC ABS, DBLE, DIMAG, MAX 174: * .. 175: * .. Statement Functions .. 176: DOUBLE PRECISION CABS1 177: * .. 178: * .. Statement Function definitions .. 179: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) 180: * .. 181: * .. Executable Statements .. 182: * 183: * Decode and test the input parameters. 184: * 185: BOTHV = LSAME( SIDE, 'B' ) 186: RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV 187: LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV 188: * 189: FROMQR = LSAME( EIGSRC, 'Q' ) 190: * 191: NOINIT = LSAME( INITV, 'N' ) 192: * 193: * Set M to the number of columns required to store the selected 194: * eigenvectors. 195: * 196: M = 0 197: DO 10 K = 1, N 198: IF( SELECT( K ) ) 199: $ M = M + 1 200: 10 CONTINUE 201: * 202: INFO = 0 203: IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN 204: INFO = -1 205: ELSE IF( .NOT.FROMQR .AND. .NOT.LSAME( EIGSRC, 'N' ) ) THEN 206: INFO = -2 207: ELSE IF( .NOT.NOINIT .AND. .NOT.LSAME( INITV, 'U' ) ) THEN 208: INFO = -3 209: ELSE IF( N.LT.0 ) THEN 210: INFO = -5 211: ELSE IF( LDH.LT.MAX( 1, N ) ) THEN 212: INFO = -7 213: ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN 214: INFO = -10 215: ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN 216: INFO = -12 217: ELSE IF( MM.LT.M ) THEN 218: INFO = -13 219: END IF 220: IF( INFO.NE.0 ) THEN 221: CALL XERBLA( 'ZHSEIN', -INFO ) 222: RETURN 223: END IF 224: * 225: * Quick return if possible. 226: * 227: IF( N.EQ.0 ) 228: $ RETURN 229: * 230: * Set machine-dependent constants. 231: * 232: UNFL = DLAMCH( 'Safe minimum' ) 233: ULP = DLAMCH( 'Precision' ) 234: SMLNUM = UNFL*( N / ULP ) 235: * 236: LDWORK = N 237: * 238: KL = 1 239: KLN = 0 240: IF( FROMQR ) THEN 241: KR = 0 242: ELSE 243: KR = N 244: END IF 245: KS = 1 246: * 247: DO 100 K = 1, N 248: IF( SELECT( K ) ) THEN 249: * 250: * Compute eigenvector(s) corresponding to W(K). 251: * 252: IF( FROMQR ) THEN 253: * 254: * If affiliation of eigenvalues is known, check whether 255: * the matrix splits. 256: * 257: * Determine KL and KR such that 1 <= KL <= K <= KR <= N 258: * and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or 259: * KR = N). 260: * 261: * Then inverse iteration can be performed with the 262: * submatrix H(KL:N,KL:N) for a left eigenvector, and with 263: * the submatrix H(1:KR,1:KR) for a right eigenvector. 264: * 265: DO 20 I = K, KL + 1, -1 266: IF( H( I, I-1 ).EQ.ZERO ) 267: $ GO TO 30 268: 20 CONTINUE 269: 30 CONTINUE 270: KL = I 271: IF( K.GT.KR ) THEN 272: DO 40 I = K, N - 1 273: IF( H( I+1, I ).EQ.ZERO ) 274: $ GO TO 50 275: 40 CONTINUE 276: 50 CONTINUE 277: KR = I 278: END IF 279: END IF 280: * 281: IF( KL.NE.KLN ) THEN 282: KLN = KL 283: * 284: * Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it 285: * has not ben computed before. 286: * 287: HNORM = ZLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, RWORK ) 288: IF( HNORM.GT.RZERO ) THEN 289: EPS3 = HNORM*ULP 290: ELSE 291: EPS3 = SMLNUM 292: END IF 293: END IF 294: * 295: * Perturb eigenvalue if it is close to any previous 296: * selected eigenvalues affiliated to the submatrix 297: * H(KL:KR,KL:KR). Close roots are modified by EPS3. 298: * 299: WK = W( K ) 300: 60 CONTINUE 301: DO 70 I = K - 1, KL, -1 302: IF( SELECT( I ) .AND. CABS1( W( I )-WK ).LT.EPS3 ) THEN 303: WK = WK + EPS3 304: GO TO 60 305: END IF 306: 70 CONTINUE 307: W( K ) = WK 308: * 309: IF( LEFTV ) THEN 310: * 311: * Compute left eigenvector. 312: * 313: CALL ZLAEIN( .FALSE., NOINIT, N-KL+1, H( KL, KL ), LDH, 314: $ WK, VL( KL, KS ), WORK, LDWORK, RWORK, EPS3, 315: $ SMLNUM, IINFO ) 316: IF( IINFO.GT.0 ) THEN 317: INFO = INFO + 1 318: IFAILL( KS ) = K 319: ELSE 320: IFAILL( KS ) = 0 321: END IF 322: DO 80 I = 1, KL - 1 323: VL( I, KS ) = ZERO 324: 80 CONTINUE 325: END IF 326: IF( RIGHTV ) THEN 327: * 328: * Compute right eigenvector. 329: * 330: CALL ZLAEIN( .TRUE., NOINIT, KR, H, LDH, WK, VR( 1, KS ), 331: $ WORK, LDWORK, RWORK, EPS3, SMLNUM, IINFO ) 332: IF( IINFO.GT.0 ) THEN 333: INFO = INFO + 1 334: IFAILR( KS ) = K 335: ELSE 336: IFAILR( KS ) = 0 337: END IF 338: DO 90 I = KR + 1, N 339: VR( I, KS ) = ZERO 340: 90 CONTINUE 341: END IF 342: KS = KS + 1 343: END IF 344: 100 CONTINUE 345: * 346: RETURN 347: * 348: * End of ZHSEIN 349: * 350: END 351: