1: *> \brief \b DLAEIN computes a specified right or left eigenvector of an upper Hessenberg matrix by inverse iteration.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLAEIN + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaein.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaein.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaein.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
22: * LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO )
23: *
24: * .. Scalar Arguments ..
25: * LOGICAL NOINIT, RIGHTV
26: * INTEGER INFO, LDB, LDH, N
27: * DOUBLE PRECISION BIGNUM, EPS3, SMLNUM, WI, WR
28: * ..
29: * .. Array Arguments ..
30: * DOUBLE PRECISION B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
31: * $ WORK( * )
32: * ..
33: *
34: *
35: *> \par Purpose:
36: * =============
37: *>
38: *> \verbatim
39: *>
40: *> DLAEIN uses inverse iteration to find a right or left eigenvector
41: *> corresponding to the eigenvalue (WR,WI) of a real upper Hessenberg
42: *> matrix H.
43: *> \endverbatim
44: *
45: * Arguments:
46: * ==========
47: *
48: *> \param[in] RIGHTV
49: *> \verbatim
50: *> RIGHTV is LOGICAL
51: *> = .TRUE. : compute right eigenvector;
52: *> = .FALSE.: compute left eigenvector.
53: *> \endverbatim
54: *>
55: *> \param[in] NOINIT
56: *> \verbatim
57: *> NOINIT is LOGICAL
58: *> = .TRUE. : no initial vector supplied in (VR,VI).
59: *> = .FALSE.: initial vector supplied in (VR,VI).
60: *> \endverbatim
61: *>
62: *> \param[in] N
63: *> \verbatim
64: *> N is INTEGER
65: *> The order of the matrix H. N >= 0.
66: *> \endverbatim
67: *>
68: *> \param[in] H
69: *> \verbatim
70: *> H is DOUBLE PRECISION array, dimension (LDH,N)
71: *> The upper Hessenberg matrix H.
72: *> \endverbatim
73: *>
74: *> \param[in] LDH
75: *> \verbatim
76: *> LDH is INTEGER
77: *> The leading dimension of the array H. LDH >= max(1,N).
78: *> \endverbatim
79: *>
80: *> \param[in] WR
81: *> \verbatim
82: *> WR is DOUBLE PRECISION
83: *> \endverbatim
84: *>
85: *> \param[in] WI
86: *> \verbatim
87: *> WI is DOUBLE PRECISION
88: *> The real and imaginary parts of the eigenvalue of H whose
89: *> corresponding right or left eigenvector is to be computed.
90: *> \endverbatim
91: *>
92: *> \param[in,out] VR
93: *> \verbatim
94: *> VR is DOUBLE PRECISION array, dimension (N)
95: *> \endverbatim
96: *>
97: *> \param[in,out] VI
98: *> \verbatim
99: *> VI is DOUBLE PRECISION array, dimension (N)
100: *> On entry, if NOINIT = .FALSE. and WI = 0.0, VR must contain
101: *> a real starting vector for inverse iteration using the real
102: *> eigenvalue WR; if NOINIT = .FALSE. and WI.ne.0.0, VR and VI
103: *> must contain the real and imaginary parts of a complex
104: *> starting vector for inverse iteration using the complex
105: *> eigenvalue (WR,WI); otherwise VR and VI need not be set.
106: *> On exit, if WI = 0.0 (real eigenvalue), VR contains the
107: *> computed real eigenvector; if WI.ne.0.0 (complex eigenvalue),
108: *> VR and VI contain the real and imaginary parts of the
109: *> computed complex eigenvector. The eigenvector is normalized
110: *> so that the component of largest magnitude has magnitude 1;
111: *> here the magnitude of a complex number (x,y) is taken to be
112: *> |x| + |y|.
113: *> VI is not referenced if WI = 0.0.
114: *> \endverbatim
115: *>
116: *> \param[out] B
117: *> \verbatim
118: *> B is DOUBLE PRECISION array, dimension (LDB,N)
119: *> \endverbatim
120: *>
121: *> \param[in] LDB
122: *> \verbatim
123: *> LDB is INTEGER
124: *> The leading dimension of the array B. LDB >= N+1.
125: *> \endverbatim
126: *>
127: *> \param[out] WORK
128: *> \verbatim
129: *> WORK is DOUBLE PRECISION array, dimension (N)
130: *> \endverbatim
131: *>
132: *> \param[in] EPS3
133: *> \verbatim
134: *> EPS3 is DOUBLE PRECISION
135: *> A small machine-dependent value which is used to perturb
136: *> close eigenvalues, and to replace zero pivots.
137: *> \endverbatim
138: *>
139: *> \param[in] SMLNUM
140: *> \verbatim
141: *> SMLNUM is DOUBLE PRECISION
142: *> A machine-dependent value close to the underflow threshold.
143: *> \endverbatim
144: *>
145: *> \param[in] BIGNUM
146: *> \verbatim
147: *> BIGNUM is DOUBLE PRECISION
148: *> A machine-dependent value close to the overflow threshold.
149: *> \endverbatim
150: *>
151: *> \param[out] INFO
152: *> \verbatim
153: *> INFO is INTEGER
154: *> = 0: successful exit
155: *> = 1: inverse iteration did not converge; VR is set to the
156: *> last iterate, and so is VI if WI.ne.0.0.
157: *> \endverbatim
158: *
159: * Authors:
160: * ========
161: *
162: *> \author Univ. of Tennessee
163: *> \author Univ. of California Berkeley
164: *> \author Univ. of Colorado Denver
165: *> \author NAG Ltd.
166: *
167: *> \ingroup doubleOTHERauxiliary
168: *
169: * =====================================================================
170: SUBROUTINE DLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
171: $ LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO )
172: *
173: * -- LAPACK auxiliary routine --
174: * -- LAPACK is a software package provided by Univ. of Tennessee, --
175: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
176: *
177: * .. Scalar Arguments ..
178: LOGICAL NOINIT, RIGHTV
179: INTEGER INFO, LDB, LDH, N
180: DOUBLE PRECISION BIGNUM, EPS3, SMLNUM, WI, WR
181: * ..
182: * .. Array Arguments ..
183: DOUBLE PRECISION B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
184: $ WORK( * )
185: * ..
186: *
187: * =====================================================================
188: *
189: * .. Parameters ..
190: DOUBLE PRECISION ZERO, ONE, TENTH
191: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TENTH = 1.0D-1 )
192: * ..
193: * .. Local Scalars ..
194: CHARACTER NORMIN, TRANS
195: INTEGER I, I1, I2, I3, IERR, ITS, J
196: DOUBLE PRECISION ABSBII, ABSBJJ, EI, EJ, GROWTO, NORM, NRMSML,
197: $ REC, ROOTN, SCALE, TEMP, VCRIT, VMAX, VNORM, W,
198: $ W1, X, XI, XR, Y
199: * ..
200: * .. External Functions ..
201: INTEGER IDAMAX
202: DOUBLE PRECISION DASUM, DLAPY2, DNRM2
203: EXTERNAL IDAMAX, DASUM, DLAPY2, DNRM2
204: * ..
205: * .. External Subroutines ..
206: EXTERNAL DLADIV, DLATRS, DSCAL
207: * ..
208: * .. Intrinsic Functions ..
209: INTRINSIC ABS, DBLE, MAX, SQRT
210: * ..
211: * .. Executable Statements ..
212: *
213: INFO = 0
214: *
215: * GROWTO is the threshold used in the acceptance test for an
216: * eigenvector.
217: *
218: ROOTN = SQRT( DBLE( N ) )
219: GROWTO = TENTH / ROOTN
220: NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM
221: *
222: * Form B = H - (WR,WI)*I (except that the subdiagonal elements and
223: * the imaginary parts of the diagonal elements are not stored).
224: *
225: DO 20 J = 1, N
226: DO 10 I = 1, J - 1
227: B( I, J ) = H( I, J )
228: 10 CONTINUE
229: B( J, J ) = H( J, J ) - WR
230: 20 CONTINUE
231: *
232: IF( WI.EQ.ZERO ) THEN
233: *
234: * Real eigenvalue.
235: *
236: IF( NOINIT ) THEN
237: *
238: * Set initial vector.
239: *
240: DO 30 I = 1, N
241: VR( I ) = EPS3
242: 30 CONTINUE
243: ELSE
244: *
245: * Scale supplied initial vector.
246: *
247: VNORM = DNRM2( N, VR, 1 )
248: CALL DSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), VR,
249: $ 1 )
250: END IF
251: *
252: IF( RIGHTV ) THEN
253: *
254: * LU decomposition with partial pivoting of B, replacing zero
255: * pivots by EPS3.
256: *
257: DO 60 I = 1, N - 1
258: EI = H( I+1, I )
259: IF( ABS( B( I, I ) ).LT.ABS( EI ) ) THEN
260: *
261: * Interchange rows and eliminate.
262: *
263: X = B( I, I ) / EI
264: B( I, I ) = EI
265: DO 40 J = I + 1, N
266: TEMP = B( I+1, J )
267: B( I+1, J ) = B( I, J ) - X*TEMP
268: B( I, J ) = TEMP
269: 40 CONTINUE
270: ELSE
271: *
272: * Eliminate without interchange.
273: *
274: IF( B( I, I ).EQ.ZERO )
275: $ B( I, I ) = EPS3
276: X = EI / B( I, I )
277: IF( X.NE.ZERO ) THEN
278: DO 50 J = I + 1, N
279: B( I+1, J ) = B( I+1, J ) - X*B( I, J )
280: 50 CONTINUE
281: END IF
282: END IF
283: 60 CONTINUE
284: IF( B( N, N ).EQ.ZERO )
285: $ B( N, N ) = EPS3
286: *
287: TRANS = 'N'
288: *
289: ELSE
290: *
291: * UL decomposition with partial pivoting of B, replacing zero
292: * pivots by EPS3.
293: *
294: DO 90 J = N, 2, -1
295: EJ = H( J, J-1 )
296: IF( ABS( B( J, J ) ).LT.ABS( EJ ) ) THEN
297: *
298: * Interchange columns and eliminate.
299: *
300: X = B( J, J ) / EJ
301: B( J, J ) = EJ
302: DO 70 I = 1, J - 1
303: TEMP = B( I, J-1 )
304: B( I, J-1 ) = B( I, J ) - X*TEMP
305: B( I, J ) = TEMP
306: 70 CONTINUE
307: ELSE
308: *
309: * Eliminate without interchange.
310: *
311: IF( B( J, J ).EQ.ZERO )
312: $ B( J, J ) = EPS3
313: X = EJ / B( J, J )
314: IF( X.NE.ZERO ) THEN
315: DO 80 I = 1, J - 1
316: B( I, J-1 ) = B( I, J-1 ) - X*B( I, J )
317: 80 CONTINUE
318: END IF
319: END IF
320: 90 CONTINUE
321: IF( B( 1, 1 ).EQ.ZERO )
322: $ B( 1, 1 ) = EPS3
323: *
324: TRANS = 'T'
325: *
326: END IF
327: *
328: NORMIN = 'N'
329: DO 110 ITS = 1, N
330: *
331: * Solve U*x = scale*v for a right eigenvector
332: * or U**T*x = scale*v for a left eigenvector,
333: * overwriting x on v.
334: *
335: CALL DLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB,
336: $ VR, SCALE, WORK, IERR )
337: NORMIN = 'Y'
338: *
339: * Test for sufficient growth in the norm of v.
340: *
341: VNORM = DASUM( N, VR, 1 )
342: IF( VNORM.GE.GROWTO*SCALE )
343: $ GO TO 120
344: *
345: * Choose new orthogonal starting vector and try again.
346: *
347: TEMP = EPS3 / ( ROOTN+ONE )
348: VR( 1 ) = EPS3
349: DO 100 I = 2, N
350: VR( I ) = TEMP
351: 100 CONTINUE
352: VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
353: 110 CONTINUE
354: *
355: * Failure to find eigenvector in N iterations.
356: *
357: INFO = 1
358: *
359: 120 CONTINUE
360: *
361: * Normalize eigenvector.
362: *
363: I = IDAMAX( N, VR, 1 )
364: CALL DSCAL( N, ONE / ABS( VR( I ) ), VR, 1 )
365: ELSE
366: *
367: * Complex eigenvalue.
368: *
369: IF( NOINIT ) THEN
370: *
371: * Set initial vector.
372: *
373: DO 130 I = 1, N
374: VR( I ) = EPS3
375: VI( I ) = ZERO
376: 130 CONTINUE
377: ELSE
378: *
379: * Scale supplied initial vector.
380: *
381: NORM = DLAPY2( DNRM2( N, VR, 1 ), DNRM2( N, VI, 1 ) )
382: REC = ( EPS3*ROOTN ) / MAX( NORM, NRMSML )
383: CALL DSCAL( N, REC, VR, 1 )
384: CALL DSCAL( N, REC, VI, 1 )
385: END IF
386: *
387: IF( RIGHTV ) THEN
388: *
389: * LU decomposition with partial pivoting of B, replacing zero
390: * pivots by EPS3.
391: *
392: * The imaginary part of the (i,j)-th element of U is stored in
393: * B(j+1,i).
394: *
395: B( 2, 1 ) = -WI
396: DO 140 I = 2, N
397: B( I+1, 1 ) = ZERO
398: 140 CONTINUE
399: *
400: DO 170 I = 1, N - 1
401: ABSBII = DLAPY2( B( I, I ), B( I+1, I ) )
402: EI = H( I+1, I )
403: IF( ABSBII.LT.ABS( EI ) ) THEN
404: *
405: * Interchange rows and eliminate.
406: *
407: XR = B( I, I ) / EI
408: XI = B( I+1, I ) / EI
409: B( I, I ) = EI
410: B( I+1, I ) = ZERO
411: DO 150 J = I + 1, N
412: TEMP = B( I+1, J )
413: B( I+1, J ) = B( I, J ) - XR*TEMP
414: B( J+1, I+1 ) = B( J+1, I ) - XI*TEMP
415: B( I, J ) = TEMP
416: B( J+1, I ) = ZERO
417: 150 CONTINUE
418: B( I+2, I ) = -WI
419: B( I+1, I+1 ) = B( I+1, I+1 ) - XI*WI
420: B( I+2, I+1 ) = B( I+2, I+1 ) + XR*WI
421: ELSE
422: *
423: * Eliminate without interchanging rows.
424: *
425: IF( ABSBII.EQ.ZERO ) THEN
426: B( I, I ) = EPS3
427: B( I+1, I ) = ZERO
428: ABSBII = EPS3
429: END IF
430: EI = ( EI / ABSBII ) / ABSBII
431: XR = B( I, I )*EI
432: XI = -B( I+1, I )*EI
433: DO 160 J = I + 1, N
434: B( I+1, J ) = B( I+1, J ) - XR*B( I, J ) +
435: $ XI*B( J+1, I )
436: B( J+1, I+1 ) = -XR*B( J+1, I ) - XI*B( I, J )
437: 160 CONTINUE
438: B( I+2, I+1 ) = B( I+2, I+1 ) - WI
439: END IF
440: *
441: * Compute 1-norm of offdiagonal elements of i-th row.
442: *
443: WORK( I ) = DASUM( N-I, B( I, I+1 ), LDB ) +
444: $ DASUM( N-I, B( I+2, I ), 1 )
445: 170 CONTINUE
446: IF( B( N, N ).EQ.ZERO .AND. B( N+1, N ).EQ.ZERO )
447: $ B( N, N ) = EPS3
448: WORK( N ) = ZERO
449: *
450: I1 = N
451: I2 = 1
452: I3 = -1
453: ELSE
454: *
455: * UL decomposition with partial pivoting of conjg(B),
456: * replacing zero pivots by EPS3.
457: *
458: * The imaginary part of the (i,j)-th element of U is stored in
459: * B(j+1,i).
460: *
461: B( N+1, N ) = WI
462: DO 180 J = 1, N - 1
463: B( N+1, J ) = ZERO
464: 180 CONTINUE
465: *
466: DO 210 J = N, 2, -1
467: EJ = H( J, J-1 )
468: ABSBJJ = DLAPY2( B( J, J ), B( J+1, J ) )
469: IF( ABSBJJ.LT.ABS( EJ ) ) THEN
470: *
471: * Interchange columns and eliminate
472: *
473: XR = B( J, J ) / EJ
474: XI = B( J+1, J ) / EJ
475: B( J, J ) = EJ
476: B( J+1, J ) = ZERO
477: DO 190 I = 1, J - 1
478: TEMP = B( I, J-1 )
479: B( I, J-1 ) = B( I, J ) - XR*TEMP
480: B( J, I ) = B( J+1, I ) - XI*TEMP
481: B( I, J ) = TEMP
482: B( J+1, I ) = ZERO
483: 190 CONTINUE
484: B( J+1, J-1 ) = WI
485: B( J-1, J-1 ) = B( J-1, J-1 ) + XI*WI
486: B( J, J-1 ) = B( J, J-1 ) - XR*WI
487: ELSE
488: *
489: * Eliminate without interchange.
490: *
491: IF( ABSBJJ.EQ.ZERO ) THEN
492: B( J, J ) = EPS3
493: B( J+1, J ) = ZERO
494: ABSBJJ = EPS3
495: END IF
496: EJ = ( EJ / ABSBJJ ) / ABSBJJ
497: XR = B( J, J )*EJ
498: XI = -B( J+1, J )*EJ
499: DO 200 I = 1, J - 1
500: B( I, J-1 ) = B( I, J-1 ) - XR*B( I, J ) +
501: $ XI*B( J+1, I )
502: B( J, I ) = -XR*B( J+1, I ) - XI*B( I, J )
503: 200 CONTINUE
504: B( J, J-1 ) = B( J, J-1 ) + WI
505: END IF
506: *
507: * Compute 1-norm of offdiagonal elements of j-th column.
508: *
509: WORK( J ) = DASUM( J-1, B( 1, J ), 1 ) +
510: $ DASUM( J-1, B( J+1, 1 ), LDB )
511: 210 CONTINUE
512: IF( B( 1, 1 ).EQ.ZERO .AND. B( 2, 1 ).EQ.ZERO )
513: $ B( 1, 1 ) = EPS3
514: WORK( 1 ) = ZERO
515: *
516: I1 = 1
517: I2 = N
518: I3 = 1
519: END IF
520: *
521: DO 270 ITS = 1, N
522: SCALE = ONE
523: VMAX = ONE
524: VCRIT = BIGNUM
525: *
526: * Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector,
527: * or U**T*(xr,xi) = scale*(vr,vi) for a left eigenvector,
528: * overwriting (xr,xi) on (vr,vi).
529: *
530: DO 250 I = I1, I2, I3
531: *
532: IF( WORK( I ).GT.VCRIT ) THEN
533: REC = ONE / VMAX
534: CALL DSCAL( N, REC, VR, 1 )
535: CALL DSCAL( N, REC, VI, 1 )
536: SCALE = SCALE*REC
537: VMAX = ONE
538: VCRIT = BIGNUM
539: END IF
540: *
541: XR = VR( I )
542: XI = VI( I )
543: IF( RIGHTV ) THEN
544: DO 220 J = I + 1, N
545: XR = XR - B( I, J )*VR( J ) + B( J+1, I )*VI( J )
546: XI = XI - B( I, J )*VI( J ) - B( J+1, I )*VR( J )
547: 220 CONTINUE
548: ELSE
549: DO 230 J = 1, I - 1
550: XR = XR - B( J, I )*VR( J ) + B( I+1, J )*VI( J )
551: XI = XI - B( J, I )*VI( J ) - B( I+1, J )*VR( J )
552: 230 CONTINUE
553: END IF
554: *
555: W = ABS( B( I, I ) ) + ABS( B( I+1, I ) )
556: IF( W.GT.SMLNUM ) THEN
557: IF( W.LT.ONE ) THEN
558: W1 = ABS( XR ) + ABS( XI )
559: IF( W1.GT.W*BIGNUM ) THEN
560: REC = ONE / W1
561: CALL DSCAL( N, REC, VR, 1 )
562: CALL DSCAL( N, REC, VI, 1 )
563: XR = VR( I )
564: XI = VI( I )
565: SCALE = SCALE*REC
566: VMAX = VMAX*REC
567: END IF
568: END IF
569: *
570: * Divide by diagonal element of B.
571: *
572: CALL DLADIV( XR, XI, B( I, I ), B( I+1, I ), VR( I ),
573: $ VI( I ) )
574: VMAX = MAX( ABS( VR( I ) )+ABS( VI( I ) ), VMAX )
575: VCRIT = BIGNUM / VMAX
576: ELSE
577: DO 240 J = 1, N
578: VR( J ) = ZERO
579: VI( J ) = ZERO
580: 240 CONTINUE
581: VR( I ) = ONE
582: VI( I ) = ONE
583: SCALE = ZERO
584: VMAX = ONE
585: VCRIT = BIGNUM
586: END IF
587: 250 CONTINUE
588: *
589: * Test for sufficient growth in the norm of (VR,VI).
590: *
591: VNORM = DASUM( N, VR, 1 ) + DASUM( N, VI, 1 )
592: IF( VNORM.GE.GROWTO*SCALE )
593: $ GO TO 280
594: *
595: * Choose a new orthogonal starting vector and try again.
596: *
597: Y = EPS3 / ( ROOTN+ONE )
598: VR( 1 ) = EPS3
599: VI( 1 ) = ZERO
600: *
601: DO 260 I = 2, N
602: VR( I ) = Y
603: VI( I ) = ZERO
604: 260 CONTINUE
605: VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
606: 270 CONTINUE
607: *
608: * Failure to find eigenvector in N iterations
609: *
610: INFO = 1
611: *
612: 280 CONTINUE
613: *
614: * Normalize eigenvector.
615: *
616: VNORM = ZERO
617: DO 290 I = 1, N
618: VNORM = MAX( VNORM, ABS( VR( I ) )+ABS( VI( I ) ) )
619: 290 CONTINUE
620: CALL DSCAL( N, ONE / VNORM, VR, 1 )
621: CALL DSCAL( N, ONE / VNORM, VI, 1 )
622: *
623: END IF
624: *
625: RETURN
626: *
627: * End of DLAEIN
628: *
629: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaein.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack