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: *> \date September 2012
168: *
169: *> \ingroup doubleOTHERauxiliary
170: *
171: * =====================================================================
172: SUBROUTINE DLAEIN( RIGHTV, NOINIT, N, H, LDH, WR, WI, VR, VI, B,
173: $ LDB, WORK, EPS3, SMLNUM, BIGNUM, INFO )
174: *
175: * -- LAPACK auxiliary routine (version 3.4.2) --
176: * -- LAPACK is a software package provided by Univ. of Tennessee, --
177: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
178: * September 2012
179: *
180: * .. Scalar Arguments ..
181: LOGICAL NOINIT, RIGHTV
182: INTEGER INFO, LDB, LDH, N
183: DOUBLE PRECISION BIGNUM, EPS3, SMLNUM, WI, WR
184: * ..
185: * .. Array Arguments ..
186: DOUBLE PRECISION B( LDB, * ), H( LDH, * ), VI( * ), VR( * ),
187: $ WORK( * )
188: * ..
189: *
190: * =====================================================================
191: *
192: * .. Parameters ..
193: DOUBLE PRECISION ZERO, ONE, TENTH
194: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TENTH = 1.0D-1 )
195: * ..
196: * .. Local Scalars ..
197: CHARACTER NORMIN, TRANS
198: INTEGER I, I1, I2, I3, IERR, ITS, J
199: DOUBLE PRECISION ABSBII, ABSBJJ, EI, EJ, GROWTO, NORM, NRMSML,
200: $ REC, ROOTN, SCALE, TEMP, VCRIT, VMAX, VNORM, W,
201: $ W1, X, XI, XR, Y
202: * ..
203: * .. External Functions ..
204: INTEGER IDAMAX
205: DOUBLE PRECISION DASUM, DLAPY2, DNRM2
206: EXTERNAL IDAMAX, DASUM, DLAPY2, DNRM2
207: * ..
208: * .. External Subroutines ..
209: EXTERNAL DLADIV, DLATRS, DSCAL
210: * ..
211: * .. Intrinsic Functions ..
212: INTRINSIC ABS, DBLE, MAX, SQRT
213: * ..
214: * .. Executable Statements ..
215: *
216: INFO = 0
217: *
218: * GROWTO is the threshold used in the acceptance test for an
219: * eigenvector.
220: *
221: ROOTN = SQRT( DBLE( N ) )
222: GROWTO = TENTH / ROOTN
223: NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM
224: *
225: * Form B = H - (WR,WI)*I (except that the subdiagonal elements and
226: * the imaginary parts of the diagonal elements are not stored).
227: *
228: DO 20 J = 1, N
229: DO 10 I = 1, J - 1
230: B( I, J ) = H( I, J )
231: 10 CONTINUE
232: B( J, J ) = H( J, J ) - WR
233: 20 CONTINUE
234: *
235: IF( WI.EQ.ZERO ) THEN
236: *
237: * Real eigenvalue.
238: *
239: IF( NOINIT ) THEN
240: *
241: * Set initial vector.
242: *
243: DO 30 I = 1, N
244: VR( I ) = EPS3
245: 30 CONTINUE
246: ELSE
247: *
248: * Scale supplied initial vector.
249: *
250: VNORM = DNRM2( N, VR, 1 )
251: CALL DSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), VR,
252: $ 1 )
253: END IF
254: *
255: IF( RIGHTV ) THEN
256: *
257: * LU decomposition with partial pivoting of B, replacing zero
258: * pivots by EPS3.
259: *
260: DO 60 I = 1, N - 1
261: EI = H( I+1, I )
262: IF( ABS( B( I, I ) ).LT.ABS( EI ) ) THEN
263: *
264: * Interchange rows and eliminate.
265: *
266: X = B( I, I ) / EI
267: B( I, I ) = EI
268: DO 40 J = I + 1, N
269: TEMP = B( I+1, J )
270: B( I+1, J ) = B( I, J ) - X*TEMP
271: B( I, J ) = TEMP
272: 40 CONTINUE
273: ELSE
274: *
275: * Eliminate without interchange.
276: *
277: IF( B( I, I ).EQ.ZERO )
278: $ B( I, I ) = EPS3
279: X = EI / B( I, I )
280: IF( X.NE.ZERO ) THEN
281: DO 50 J = I + 1, N
282: B( I+1, J ) = B( I+1, J ) - X*B( I, J )
283: 50 CONTINUE
284: END IF
285: END IF
286: 60 CONTINUE
287: IF( B( N, N ).EQ.ZERO )
288: $ B( N, N ) = EPS3
289: *
290: TRANS = 'N'
291: *
292: ELSE
293: *
294: * UL decomposition with partial pivoting of B, replacing zero
295: * pivots by EPS3.
296: *
297: DO 90 J = N, 2, -1
298: EJ = H( J, J-1 )
299: IF( ABS( B( J, J ) ).LT.ABS( EJ ) ) THEN
300: *
301: * Interchange columns and eliminate.
302: *
303: X = B( J, J ) / EJ
304: B( J, J ) = EJ
305: DO 70 I = 1, J - 1
306: TEMP = B( I, J-1 )
307: B( I, J-1 ) = B( I, J ) - X*TEMP
308: B( I, J ) = TEMP
309: 70 CONTINUE
310: ELSE
311: *
312: * Eliminate without interchange.
313: *
314: IF( B( J, J ).EQ.ZERO )
315: $ B( J, J ) = EPS3
316: X = EJ / B( J, J )
317: IF( X.NE.ZERO ) THEN
318: DO 80 I = 1, J - 1
319: B( I, J-1 ) = B( I, J-1 ) - X*B( I, J )
320: 80 CONTINUE
321: END IF
322: END IF
323: 90 CONTINUE
324: IF( B( 1, 1 ).EQ.ZERO )
325: $ B( 1, 1 ) = EPS3
326: *
327: TRANS = 'T'
328: *
329: END IF
330: *
331: NORMIN = 'N'
332: DO 110 ITS = 1, N
333: *
334: * Solve U*x = scale*v for a right eigenvector
335: * or U**T*x = scale*v for a left eigenvector,
336: * overwriting x on v.
337: *
338: CALL DLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB,
339: $ VR, SCALE, WORK, IERR )
340: NORMIN = 'Y'
341: *
342: * Test for sufficient growth in the norm of v.
343: *
344: VNORM = DASUM( N, VR, 1 )
345: IF( VNORM.GE.GROWTO*SCALE )
346: $ GO TO 120
347: *
348: * Choose new orthogonal starting vector and try again.
349: *
350: TEMP = EPS3 / ( ROOTN+ONE )
351: VR( 1 ) = EPS3
352: DO 100 I = 2, N
353: VR( I ) = TEMP
354: 100 CONTINUE
355: VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
356: 110 CONTINUE
357: *
358: * Failure to find eigenvector in N iterations.
359: *
360: INFO = 1
361: *
362: 120 CONTINUE
363: *
364: * Normalize eigenvector.
365: *
366: I = IDAMAX( N, VR, 1 )
367: CALL DSCAL( N, ONE / ABS( VR( I ) ), VR, 1 )
368: ELSE
369: *
370: * Complex eigenvalue.
371: *
372: IF( NOINIT ) THEN
373: *
374: * Set initial vector.
375: *
376: DO 130 I = 1, N
377: VR( I ) = EPS3
378: VI( I ) = ZERO
379: 130 CONTINUE
380: ELSE
381: *
382: * Scale supplied initial vector.
383: *
384: NORM = DLAPY2( DNRM2( N, VR, 1 ), DNRM2( N, VI, 1 ) )
385: REC = ( EPS3*ROOTN ) / MAX( NORM, NRMSML )
386: CALL DSCAL( N, REC, VR, 1 )
387: CALL DSCAL( N, REC, VI, 1 )
388: END IF
389: *
390: IF( RIGHTV ) THEN
391: *
392: * LU decomposition with partial pivoting of B, replacing zero
393: * pivots by EPS3.
394: *
395: * The imaginary part of the (i,j)-th element of U is stored in
396: * B(j+1,i).
397: *
398: B( 2, 1 ) = -WI
399: DO 140 I = 2, N
400: B( I+1, 1 ) = ZERO
401: 140 CONTINUE
402: *
403: DO 170 I = 1, N - 1
404: ABSBII = DLAPY2( B( I, I ), B( I+1, I ) )
405: EI = H( I+1, I )
406: IF( ABSBII.LT.ABS( EI ) ) THEN
407: *
408: * Interchange rows and eliminate.
409: *
410: XR = B( I, I ) / EI
411: XI = B( I+1, I ) / EI
412: B( I, I ) = EI
413: B( I+1, I ) = ZERO
414: DO 150 J = I + 1, N
415: TEMP = B( I+1, J )
416: B( I+1, J ) = B( I, J ) - XR*TEMP
417: B( J+1, I+1 ) = B( J+1, I ) - XI*TEMP
418: B( I, J ) = TEMP
419: B( J+1, I ) = ZERO
420: 150 CONTINUE
421: B( I+2, I ) = -WI
422: B( I+1, I+1 ) = B( I+1, I+1 ) - XI*WI
423: B( I+2, I+1 ) = B( I+2, I+1 ) + XR*WI
424: ELSE
425: *
426: * Eliminate without interchanging rows.
427: *
428: IF( ABSBII.EQ.ZERO ) THEN
429: B( I, I ) = EPS3
430: B( I+1, I ) = ZERO
431: ABSBII = EPS3
432: END IF
433: EI = ( EI / ABSBII ) / ABSBII
434: XR = B( I, I )*EI
435: XI = -B( I+1, I )*EI
436: DO 160 J = I + 1, N
437: B( I+1, J ) = B( I+1, J ) - XR*B( I, J ) +
438: $ XI*B( J+1, I )
439: B( J+1, I+1 ) = -XR*B( J+1, I ) - XI*B( I, J )
440: 160 CONTINUE
441: B( I+2, I+1 ) = B( I+2, I+1 ) - WI
442: END IF
443: *
444: * Compute 1-norm of offdiagonal elements of i-th row.
445: *
446: WORK( I ) = DASUM( N-I, B( I, I+1 ), LDB ) +
447: $ DASUM( N-I, B( I+2, I ), 1 )
448: 170 CONTINUE
449: IF( B( N, N ).EQ.ZERO .AND. B( N+1, N ).EQ.ZERO )
450: $ B( N, N ) = EPS3
451: WORK( N ) = ZERO
452: *
453: I1 = N
454: I2 = 1
455: I3 = -1
456: ELSE
457: *
458: * UL decomposition with partial pivoting of conjg(B),
459: * replacing zero pivots by EPS3.
460: *
461: * The imaginary part of the (i,j)-th element of U is stored in
462: * B(j+1,i).
463: *
464: B( N+1, N ) = WI
465: DO 180 J = 1, N - 1
466: B( N+1, J ) = ZERO
467: 180 CONTINUE
468: *
469: DO 210 J = N, 2, -1
470: EJ = H( J, J-1 )
471: ABSBJJ = DLAPY2( B( J, J ), B( J+1, J ) )
472: IF( ABSBJJ.LT.ABS( EJ ) ) THEN
473: *
474: * Interchange columns and eliminate
475: *
476: XR = B( J, J ) / EJ
477: XI = B( J+1, J ) / EJ
478: B( J, J ) = EJ
479: B( J+1, J ) = ZERO
480: DO 190 I = 1, J - 1
481: TEMP = B( I, J-1 )
482: B( I, J-1 ) = B( I, J ) - XR*TEMP
483: B( J, I ) = B( J+1, I ) - XI*TEMP
484: B( I, J ) = TEMP
485: B( J+1, I ) = ZERO
486: 190 CONTINUE
487: B( J+1, J-1 ) = WI
488: B( J-1, J-1 ) = B( J-1, J-1 ) + XI*WI
489: B( J, J-1 ) = B( J, J-1 ) - XR*WI
490: ELSE
491: *
492: * Eliminate without interchange.
493: *
494: IF( ABSBJJ.EQ.ZERO ) THEN
495: B( J, J ) = EPS3
496: B( J+1, J ) = ZERO
497: ABSBJJ = EPS3
498: END IF
499: EJ = ( EJ / ABSBJJ ) / ABSBJJ
500: XR = B( J, J )*EJ
501: XI = -B( J+1, J )*EJ
502: DO 200 I = 1, J - 1
503: B( I, J-1 ) = B( I, J-1 ) - XR*B( I, J ) +
504: $ XI*B( J+1, I )
505: B( J, I ) = -XR*B( J+1, I ) - XI*B( I, J )
506: 200 CONTINUE
507: B( J, J-1 ) = B( J, J-1 ) + WI
508: END IF
509: *
510: * Compute 1-norm of offdiagonal elements of j-th column.
511: *
512: WORK( J ) = DASUM( J-1, B( 1, J ), 1 ) +
513: $ DASUM( J-1, B( J+1, 1 ), LDB )
514: 210 CONTINUE
515: IF( B( 1, 1 ).EQ.ZERO .AND. B( 2, 1 ).EQ.ZERO )
516: $ B( 1, 1 ) = EPS3
517: WORK( 1 ) = ZERO
518: *
519: I1 = 1
520: I2 = N
521: I3 = 1
522: END IF
523: *
524: DO 270 ITS = 1, N
525: SCALE = ONE
526: VMAX = ONE
527: VCRIT = BIGNUM
528: *
529: * Solve U*(xr,xi) = scale*(vr,vi) for a right eigenvector,
530: * or U**T*(xr,xi) = scale*(vr,vi) for a left eigenvector,
531: * overwriting (xr,xi) on (vr,vi).
532: *
533: DO 250 I = I1, I2, I3
534: *
535: IF( WORK( I ).GT.VCRIT ) THEN
536: REC = ONE / VMAX
537: CALL DSCAL( N, REC, VR, 1 )
538: CALL DSCAL( N, REC, VI, 1 )
539: SCALE = SCALE*REC
540: VMAX = ONE
541: VCRIT = BIGNUM
542: END IF
543: *
544: XR = VR( I )
545: XI = VI( I )
546: IF( RIGHTV ) THEN
547: DO 220 J = I + 1, N
548: XR = XR - B( I, J )*VR( J ) + B( J+1, I )*VI( J )
549: XI = XI - B( I, J )*VI( J ) - B( J+1, I )*VR( J )
550: 220 CONTINUE
551: ELSE
552: DO 230 J = 1, I - 1
553: XR = XR - B( J, I )*VR( J ) + B( I+1, J )*VI( J )
554: XI = XI - B( J, I )*VI( J ) - B( I+1, J )*VR( J )
555: 230 CONTINUE
556: END IF
557: *
558: W = ABS( B( I, I ) ) + ABS( B( I+1, I ) )
559: IF( W.GT.SMLNUM ) THEN
560: IF( W.LT.ONE ) THEN
561: W1 = ABS( XR ) + ABS( XI )
562: IF( W1.GT.W*BIGNUM ) THEN
563: REC = ONE / W1
564: CALL DSCAL( N, REC, VR, 1 )
565: CALL DSCAL( N, REC, VI, 1 )
566: XR = VR( I )
567: XI = VI( I )
568: SCALE = SCALE*REC
569: VMAX = VMAX*REC
570: END IF
571: END IF
572: *
573: * Divide by diagonal element of B.
574: *
575: CALL DLADIV( XR, XI, B( I, I ), B( I+1, I ), VR( I ),
576: $ VI( I ) )
577: VMAX = MAX( ABS( VR( I ) )+ABS( VI( I ) ), VMAX )
578: VCRIT = BIGNUM / VMAX
579: ELSE
580: DO 240 J = 1, N
581: VR( J ) = ZERO
582: VI( J ) = ZERO
583: 240 CONTINUE
584: VR( I ) = ONE
585: VI( I ) = ONE
586: SCALE = ZERO
587: VMAX = ONE
588: VCRIT = BIGNUM
589: END IF
590: 250 CONTINUE
591: *
592: * Test for sufficient growth in the norm of (VR,VI).
593: *
594: VNORM = DASUM( N, VR, 1 ) + DASUM( N, VI, 1 )
595: IF( VNORM.GE.GROWTO*SCALE )
596: $ GO TO 280
597: *
598: * Choose a new orthogonal starting vector and try again.
599: *
600: Y = EPS3 / ( ROOTN+ONE )
601: VR( 1 ) = EPS3
602: VI( 1 ) = ZERO
603: *
604: DO 260 I = 2, N
605: VR( I ) = Y
606: VI( I ) = ZERO
607: 260 CONTINUE
608: VR( N-ITS+1 ) = VR( N-ITS+1 ) - EPS3*ROOTN
609: 270 CONTINUE
610: *
611: * Failure to find eigenvector in N iterations
612: *
613: INFO = 1
614: *
615: 280 CONTINUE
616: *
617: * Normalize eigenvector.
618: *
619: VNORM = ZERO
620: DO 290 I = 1, N
621: VNORM = MAX( VNORM, ABS( VR( I ) )+ABS( VI( I ) ) )
622: 290 CONTINUE
623: CALL DSCAL( N, ONE / VNORM, VR, 1 )
624: CALL DSCAL( N, ONE / VNORM, VI, 1 )
625: *
626: END IF
627: *
628: RETURN
629: *
630: * End of DLAEIN
631: *
632: 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