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