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