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