Annotation of rpl/lapack/lapack/zhsein.f, revision 1.8
1.8 ! bertrand 1: *> \brief \b ZHSEIN
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZHSEIN + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhsein.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhsein.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhsein.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL,
! 22: * LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL,
! 23: * IFAILR, INFO )
! 24: *
! 25: * .. Scalar Arguments ..
! 26: * CHARACTER EIGSRC, INITV, SIDE
! 27: * INTEGER INFO, LDH, LDVL, LDVR, M, MM, N
! 28: * ..
! 29: * .. Array Arguments ..
! 30: * LOGICAL SELECT( * )
! 31: * INTEGER IFAILL( * ), IFAILR( * )
! 32: * DOUBLE PRECISION RWORK( * )
! 33: * COMPLEX*16 H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
! 34: * $ W( * ), WORK( * )
! 35: * ..
! 36: *
! 37: *
! 38: *> \par Purpose:
! 39: * =============
! 40: *>
! 41: *> \verbatim
! 42: *>
! 43: *> ZHSEIN uses inverse iteration to find specified right and/or left
! 44: *> eigenvectors of a complex upper Hessenberg matrix H.
! 45: *>
! 46: *> The right eigenvector x and the left eigenvector y of the matrix H
! 47: *> corresponding to an eigenvalue w are defined by:
! 48: *>
! 49: *> H * x = w * x, y**h * H = w * y**h
! 50: *>
! 51: *> where y**h denotes the conjugate transpose of the vector y.
! 52: *> \endverbatim
! 53: *
! 54: * Arguments:
! 55: * ==========
! 56: *
! 57: *> \param[in] SIDE
! 58: *> \verbatim
! 59: *> SIDE is CHARACTER*1
! 60: *> = 'R': compute right eigenvectors only;
! 61: *> = 'L': compute left eigenvectors only;
! 62: *> = 'B': compute both right and left eigenvectors.
! 63: *> \endverbatim
! 64: *>
! 65: *> \param[in] EIGSRC
! 66: *> \verbatim
! 67: *> EIGSRC is CHARACTER*1
! 68: *> Specifies the source of eigenvalues supplied in W:
! 69: *> = 'Q': the eigenvalues were found using ZHSEQR; thus, if
! 70: *> H has zero subdiagonal elements, and so is
! 71: *> block-triangular, then the j-th eigenvalue can be
! 72: *> assumed to be an eigenvalue of the block containing
! 73: *> the j-th row/column. This property allows ZHSEIN to
! 74: *> perform inverse iteration on just one diagonal block.
! 75: *> = 'N': no assumptions are made on the correspondence
! 76: *> between eigenvalues and diagonal blocks. In this
! 77: *> case, ZHSEIN must always perform inverse iteration
! 78: *> using the whole matrix H.
! 79: *> \endverbatim
! 80: *>
! 81: *> \param[in] INITV
! 82: *> \verbatim
! 83: *> INITV is CHARACTER*1
! 84: *> = 'N': no initial vectors are supplied;
! 85: *> = 'U': user-supplied initial vectors are stored in the arrays
! 86: *> VL and/or VR.
! 87: *> \endverbatim
! 88: *>
! 89: *> \param[in] SELECT
! 90: *> \verbatim
! 91: *> SELECT is LOGICAL array, dimension (N)
! 92: *> Specifies the eigenvectors to be computed. To select the
! 93: *> eigenvector corresponding to the eigenvalue W(j),
! 94: *> SELECT(j) must be set to .TRUE..
! 95: *> \endverbatim
! 96: *>
! 97: *> \param[in] N
! 98: *> \verbatim
! 99: *> N is INTEGER
! 100: *> The order of the matrix H. N >= 0.
! 101: *> \endverbatim
! 102: *>
! 103: *> \param[in] H
! 104: *> \verbatim
! 105: *> H is COMPLEX*16 array, dimension (LDH,N)
! 106: *> The upper Hessenberg matrix H.
! 107: *> \endverbatim
! 108: *>
! 109: *> \param[in] LDH
! 110: *> \verbatim
! 111: *> LDH is INTEGER
! 112: *> The leading dimension of the array H. LDH >= max(1,N).
! 113: *> \endverbatim
! 114: *>
! 115: *> \param[in,out] W
! 116: *> \verbatim
! 117: *> W is COMPLEX*16 array, dimension (N)
! 118: *> On entry, the eigenvalues of H.
! 119: *> On exit, the real parts of W may have been altered since
! 120: *> close eigenvalues are perturbed slightly in searching for
! 121: *> independent eigenvectors.
! 122: *> \endverbatim
! 123: *>
! 124: *> \param[in,out] VL
! 125: *> \verbatim
! 126: *> VL is COMPLEX*16 array, dimension (LDVL,MM)
! 127: *> On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must
! 128: *> contain starting vectors for the inverse iteration for the
! 129: *> left eigenvectors; the starting vector for each eigenvector
! 130: *> must be in the same column in which the eigenvector will be
! 131: *> stored.
! 132: *> On exit, if SIDE = 'L' or 'B', the left eigenvectors
! 133: *> specified by SELECT will be stored consecutively in the
! 134: *> columns of VL, in the same order as their eigenvalues.
! 135: *> If SIDE = 'R', VL is not referenced.
! 136: *> \endverbatim
! 137: *>
! 138: *> \param[in] LDVL
! 139: *> \verbatim
! 140: *> LDVL is INTEGER
! 141: *> The leading dimension of the array VL.
! 142: *> LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
! 143: *> \endverbatim
! 144: *>
! 145: *> \param[in,out] VR
! 146: *> \verbatim
! 147: *> VR is COMPLEX*16 array, dimension (LDVR,MM)
! 148: *> On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must
! 149: *> contain starting vectors for the inverse iteration for the
! 150: *> right eigenvectors; the starting vector for each eigenvector
! 151: *> must be in the same column in which the eigenvector will be
! 152: *> stored.
! 153: *> On exit, if SIDE = 'R' or 'B', the right eigenvectors
! 154: *> specified by SELECT will be stored consecutively in the
! 155: *> columns of VR, in the same order as their eigenvalues.
! 156: *> If SIDE = 'L', VR is not referenced.
! 157: *> \endverbatim
! 158: *>
! 159: *> \param[in] LDVR
! 160: *> \verbatim
! 161: *> LDVR is INTEGER
! 162: *> The leading dimension of the array VR.
! 163: *> LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
! 164: *> \endverbatim
! 165: *>
! 166: *> \param[in] MM
! 167: *> \verbatim
! 168: *> MM is INTEGER
! 169: *> The number of columns in the arrays VL and/or VR. MM >= M.
! 170: *> \endverbatim
! 171: *>
! 172: *> \param[out] M
! 173: *> \verbatim
! 174: *> M is INTEGER
! 175: *> The number of columns in the arrays VL and/or VR required to
! 176: *> store the eigenvectors (= the number of .TRUE. elements in
! 177: *> SELECT).
! 178: *> \endverbatim
! 179: *>
! 180: *> \param[out] WORK
! 181: *> \verbatim
! 182: *> WORK is COMPLEX*16 array, dimension (N*N)
! 183: *> \endverbatim
! 184: *>
! 185: *> \param[out] RWORK
! 186: *> \verbatim
! 187: *> RWORK is DOUBLE PRECISION array, dimension (N)
! 188: *> \endverbatim
! 189: *>
! 190: *> \param[out] IFAILL
! 191: *> \verbatim
! 192: *> IFAILL is INTEGER array, dimension (MM)
! 193: *> If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left
! 194: *> eigenvector in the i-th column of VL (corresponding to the
! 195: *> eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the
! 196: *> eigenvector converged satisfactorily.
! 197: *> If SIDE = 'R', IFAILL is not referenced.
! 198: *> \endverbatim
! 199: *>
! 200: *> \param[out] IFAILR
! 201: *> \verbatim
! 202: *> IFAILR is INTEGER array, dimension (MM)
! 203: *> If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right
! 204: *> eigenvector in the i-th column of VR (corresponding to the
! 205: *> eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the
! 206: *> eigenvector converged satisfactorily.
! 207: *> If SIDE = 'L', IFAILR is not referenced.
! 208: *> \endverbatim
! 209: *>
! 210: *> \param[out] INFO
! 211: *> \verbatim
! 212: *> INFO is INTEGER
! 213: *> = 0: successful exit
! 214: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 215: *> > 0: if INFO = i, i is the number of eigenvectors which
! 216: *> failed to converge; see IFAILL and IFAILR for further
! 217: *> details.
! 218: *> \endverbatim
! 219: *
! 220: * Authors:
! 221: * ========
! 222: *
! 223: *> \author Univ. of Tennessee
! 224: *> \author Univ. of California Berkeley
! 225: *> \author Univ. of Colorado Denver
! 226: *> \author NAG Ltd.
! 227: *
! 228: *> \date November 2011
! 229: *
! 230: *> \ingroup complex16OTHERcomputational
! 231: *
! 232: *> \par Further Details:
! 233: * =====================
! 234: *>
! 235: *> \verbatim
! 236: *>
! 237: *> Each eigenvector is normalized so that the element of largest
! 238: *> magnitude has magnitude 1; here the magnitude of a complex number
! 239: *> (x,y) is taken to be |x|+|y|.
! 240: *> \endverbatim
! 241: *>
! 242: * =====================================================================
1.1 bertrand 243: SUBROUTINE ZHSEIN( SIDE, EIGSRC, INITV, SELECT, N, H, LDH, W, VL,
244: $ LDVL, VR, LDVR, MM, M, WORK, RWORK, IFAILL,
245: $ IFAILR, INFO )
246: *
1.8 ! bertrand 247: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 248: * -- LAPACK is a software package provided by Univ. of Tennessee, --
249: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 250: * November 2011
1.1 bertrand 251: *
252: * .. Scalar Arguments ..
253: CHARACTER EIGSRC, INITV, SIDE
254: INTEGER INFO, LDH, LDVL, LDVR, M, MM, N
255: * ..
256: * .. Array Arguments ..
257: LOGICAL SELECT( * )
258: INTEGER IFAILL( * ), IFAILR( * )
259: DOUBLE PRECISION RWORK( * )
260: COMPLEX*16 H( LDH, * ), VL( LDVL, * ), VR( LDVR, * ),
261: $ W( * ), WORK( * )
262: * ..
263: *
264: * =====================================================================
265: *
266: * .. Parameters ..
267: COMPLEX*16 ZERO
268: PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
269: DOUBLE PRECISION RZERO
270: PARAMETER ( RZERO = 0.0D+0 )
271: * ..
272: * .. Local Scalars ..
273: LOGICAL BOTHV, FROMQR, LEFTV, NOINIT, RIGHTV
274: INTEGER I, IINFO, K, KL, KLN, KR, KS, LDWORK
275: DOUBLE PRECISION EPS3, HNORM, SMLNUM, ULP, UNFL
276: COMPLEX*16 CDUM, WK
277: * ..
278: * .. External Functions ..
279: LOGICAL LSAME
280: DOUBLE PRECISION DLAMCH, ZLANHS
281: EXTERNAL LSAME, DLAMCH, ZLANHS
282: * ..
283: * .. External Subroutines ..
284: EXTERNAL XERBLA, ZLAEIN
285: * ..
286: * .. Intrinsic Functions ..
287: INTRINSIC ABS, DBLE, DIMAG, MAX
288: * ..
289: * .. Statement Functions ..
290: DOUBLE PRECISION CABS1
291: * ..
292: * .. Statement Function definitions ..
293: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
294: * ..
295: * .. Executable Statements ..
296: *
297: * Decode and test the input parameters.
298: *
299: BOTHV = LSAME( SIDE, 'B' )
300: RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
301: LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
302: *
303: FROMQR = LSAME( EIGSRC, 'Q' )
304: *
305: NOINIT = LSAME( INITV, 'N' )
306: *
307: * Set M to the number of columns required to store the selected
308: * eigenvectors.
309: *
310: M = 0
311: DO 10 K = 1, N
312: IF( SELECT( K ) )
313: $ M = M + 1
314: 10 CONTINUE
315: *
316: INFO = 0
317: IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
318: INFO = -1
319: ELSE IF( .NOT.FROMQR .AND. .NOT.LSAME( EIGSRC, 'N' ) ) THEN
320: INFO = -2
321: ELSE IF( .NOT.NOINIT .AND. .NOT.LSAME( INITV, 'U' ) ) THEN
322: INFO = -3
323: ELSE IF( N.LT.0 ) THEN
324: INFO = -5
325: ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
326: INFO = -7
327: ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
328: INFO = -10
329: ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
330: INFO = -12
331: ELSE IF( MM.LT.M ) THEN
332: INFO = -13
333: END IF
334: IF( INFO.NE.0 ) THEN
335: CALL XERBLA( 'ZHSEIN', -INFO )
336: RETURN
337: END IF
338: *
339: * Quick return if possible.
340: *
341: IF( N.EQ.0 )
342: $ RETURN
343: *
344: * Set machine-dependent constants.
345: *
346: UNFL = DLAMCH( 'Safe minimum' )
347: ULP = DLAMCH( 'Precision' )
348: SMLNUM = UNFL*( N / ULP )
349: *
350: LDWORK = N
351: *
352: KL = 1
353: KLN = 0
354: IF( FROMQR ) THEN
355: KR = 0
356: ELSE
357: KR = N
358: END IF
359: KS = 1
360: *
361: DO 100 K = 1, N
362: IF( SELECT( K ) ) THEN
363: *
364: * Compute eigenvector(s) corresponding to W(K).
365: *
366: IF( FROMQR ) THEN
367: *
368: * If affiliation of eigenvalues is known, check whether
369: * the matrix splits.
370: *
371: * Determine KL and KR such that 1 <= KL <= K <= KR <= N
372: * and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or
373: * KR = N).
374: *
375: * Then inverse iteration can be performed with the
376: * submatrix H(KL:N,KL:N) for a left eigenvector, and with
377: * the submatrix H(1:KR,1:KR) for a right eigenvector.
378: *
379: DO 20 I = K, KL + 1, -1
380: IF( H( I, I-1 ).EQ.ZERO )
381: $ GO TO 30
382: 20 CONTINUE
383: 30 CONTINUE
384: KL = I
385: IF( K.GT.KR ) THEN
386: DO 40 I = K, N - 1
387: IF( H( I+1, I ).EQ.ZERO )
388: $ GO TO 50
389: 40 CONTINUE
390: 50 CONTINUE
391: KR = I
392: END IF
393: END IF
394: *
395: IF( KL.NE.KLN ) THEN
396: KLN = KL
397: *
398: * Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it
399: * has not ben computed before.
400: *
401: HNORM = ZLANHS( 'I', KR-KL+1, H( KL, KL ), LDH, RWORK )
402: IF( HNORM.GT.RZERO ) THEN
403: EPS3 = HNORM*ULP
404: ELSE
405: EPS3 = SMLNUM
406: END IF
407: END IF
408: *
409: * Perturb eigenvalue if it is close to any previous
410: * selected eigenvalues affiliated to the submatrix
411: * H(KL:KR,KL:KR). Close roots are modified by EPS3.
412: *
413: WK = W( K )
414: 60 CONTINUE
415: DO 70 I = K - 1, KL, -1
416: IF( SELECT( I ) .AND. CABS1( W( I )-WK ).LT.EPS3 ) THEN
417: WK = WK + EPS3
418: GO TO 60
419: END IF
420: 70 CONTINUE
421: W( K ) = WK
422: *
423: IF( LEFTV ) THEN
424: *
425: * Compute left eigenvector.
426: *
427: CALL ZLAEIN( .FALSE., NOINIT, N-KL+1, H( KL, KL ), LDH,
428: $ WK, VL( KL, KS ), WORK, LDWORK, RWORK, EPS3,
429: $ SMLNUM, IINFO )
430: IF( IINFO.GT.0 ) THEN
431: INFO = INFO + 1
432: IFAILL( KS ) = K
433: ELSE
434: IFAILL( KS ) = 0
435: END IF
436: DO 80 I = 1, KL - 1
437: VL( I, KS ) = ZERO
438: 80 CONTINUE
439: END IF
440: IF( RIGHTV ) THEN
441: *
442: * Compute right eigenvector.
443: *
444: CALL ZLAEIN( .TRUE., NOINIT, KR, H, LDH, WK, VR( 1, KS ),
445: $ WORK, LDWORK, RWORK, EPS3, SMLNUM, IINFO )
446: IF( IINFO.GT.0 ) THEN
447: INFO = INFO + 1
448: IFAILR( KS ) = K
449: ELSE
450: IFAILR( KS ) = 0
451: END IF
452: DO 90 I = KR + 1, N
453: VR( I, KS ) = ZERO
454: 90 CONTINUE
455: END IF
456: KS = KS + 1
457: END IF
458: 100 CONTINUE
459: *
460: RETURN
461: *
462: * End of ZHSEIN
463: *
464: 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