Annotation of rpl/lapack/lapack/zlaein.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b ZLAEIN
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZLAEIN + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaein.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaein.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaein.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK,
! 22: * EPS3, SMLNUM, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * LOGICAL NOINIT, RIGHTV
! 26: * INTEGER INFO, LDB, LDH, N
! 27: * DOUBLE PRECISION EPS3, SMLNUM
! 28: * COMPLEX*16 W
! 29: * ..
! 30: * .. Array Arguments ..
! 31: * DOUBLE PRECISION RWORK( * )
! 32: * COMPLEX*16 B( LDB, * ), H( LDH, * ), V( * )
! 33: * ..
! 34: *
! 35: *
! 36: *> \par Purpose:
! 37: * =============
! 38: *>
! 39: *> \verbatim
! 40: *>
! 41: *> ZLAEIN uses inverse iteration to find a right or left eigenvector
! 42: *> corresponding to the eigenvalue W of a complex upper Hessenberg
! 43: *> matrix H.
! 44: *> \endverbatim
! 45: *
! 46: * Arguments:
! 47: * ==========
! 48: *
! 49: *> \param[in] RIGHTV
! 50: *> \verbatim
! 51: *> RIGHTV is LOGICAL
! 52: *> = .TRUE. : compute right eigenvector;
! 53: *> = .FALSE.: compute left eigenvector.
! 54: *> \endverbatim
! 55: *>
! 56: *> \param[in] NOINIT
! 57: *> \verbatim
! 58: *> NOINIT is LOGICAL
! 59: *> = .TRUE. : no initial vector supplied in V
! 60: *> = .FALSE.: initial vector supplied in V.
! 61: *> \endverbatim
! 62: *>
! 63: *> \param[in] N
! 64: *> \verbatim
! 65: *> N is INTEGER
! 66: *> The order of the matrix H. N >= 0.
! 67: *> \endverbatim
! 68: *>
! 69: *> \param[in] H
! 70: *> \verbatim
! 71: *> H is COMPLEX*16 array, dimension (LDH,N)
! 72: *> The upper Hessenberg matrix H.
! 73: *> \endverbatim
! 74: *>
! 75: *> \param[in] LDH
! 76: *> \verbatim
! 77: *> LDH is INTEGER
! 78: *> The leading dimension of the array H. LDH >= max(1,N).
! 79: *> \endverbatim
! 80: *>
! 81: *> \param[in] W
! 82: *> \verbatim
! 83: *> W is COMPLEX*16
! 84: *> The eigenvalue of H whose corresponding right or left
! 85: *> eigenvector is to be computed.
! 86: *> \endverbatim
! 87: *>
! 88: *> \param[in,out] V
! 89: *> \verbatim
! 90: *> V is COMPLEX*16 array, dimension (N)
! 91: *> On entry, if NOINIT = .FALSE., V must contain a starting
! 92: *> vector for inverse iteration; otherwise V need not be set.
! 93: *> On exit, V contains the computed eigenvector, normalized so
! 94: *> that the component of largest magnitude has magnitude 1; here
! 95: *> the magnitude of a complex number (x,y) is taken to be
! 96: *> |x| + |y|.
! 97: *> \endverbatim
! 98: *>
! 99: *> \param[out] B
! 100: *> \verbatim
! 101: *> B is COMPLEX*16 array, dimension (LDB,N)
! 102: *> \endverbatim
! 103: *>
! 104: *> \param[in] LDB
! 105: *> \verbatim
! 106: *> LDB is INTEGER
! 107: *> The leading dimension of the array B. LDB >= max(1,N).
! 108: *> \endverbatim
! 109: *>
! 110: *> \param[out] RWORK
! 111: *> \verbatim
! 112: *> RWORK is DOUBLE PRECISION array, dimension (N)
! 113: *> \endverbatim
! 114: *>
! 115: *> \param[in] EPS3
! 116: *> \verbatim
! 117: *> EPS3 is DOUBLE PRECISION
! 118: *> A small machine-dependent value which is used to perturb
! 119: *> close eigenvalues, and to replace zero pivots.
! 120: *> \endverbatim
! 121: *>
! 122: *> \param[in] SMLNUM
! 123: *> \verbatim
! 124: *> SMLNUM is DOUBLE PRECISION
! 125: *> A machine-dependent value close to the underflow threshold.
! 126: *> \endverbatim
! 127: *>
! 128: *> \param[out] INFO
! 129: *> \verbatim
! 130: *> INFO is INTEGER
! 131: *> = 0: successful exit
! 132: *> = 1: inverse iteration did not converge; V is set to the
! 133: *> last iterate.
! 134: *> \endverbatim
! 135: *
! 136: * Authors:
! 137: * ========
! 138: *
! 139: *> \author Univ. of Tennessee
! 140: *> \author Univ. of California Berkeley
! 141: *> \author Univ. of Colorado Denver
! 142: *> \author NAG Ltd.
! 143: *
! 144: *> \date November 2011
! 145: *
! 146: *> \ingroup complex16OTHERauxiliary
! 147: *
! 148: * =====================================================================
1.1 bertrand 149: SUBROUTINE ZLAEIN( RIGHTV, NOINIT, N, H, LDH, W, V, B, LDB, RWORK,
150: $ EPS3, SMLNUM, INFO )
151: *
1.9 ! bertrand 152: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 153: * -- LAPACK is a software package provided by Univ. of Tennessee, --
154: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 155: * November 2011
1.1 bertrand 156: *
157: * .. Scalar Arguments ..
158: LOGICAL NOINIT, RIGHTV
159: INTEGER INFO, LDB, LDH, N
160: DOUBLE PRECISION EPS3, SMLNUM
161: COMPLEX*16 W
162: * ..
163: * .. Array Arguments ..
164: DOUBLE PRECISION RWORK( * )
165: COMPLEX*16 B( LDB, * ), H( LDH, * ), V( * )
166: * ..
167: *
168: * =====================================================================
169: *
170: * .. Parameters ..
171: DOUBLE PRECISION ONE, TENTH
172: PARAMETER ( ONE = 1.0D+0, TENTH = 1.0D-1 )
173: COMPLEX*16 ZERO
174: PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
175: * ..
176: * .. Local Scalars ..
177: CHARACTER NORMIN, TRANS
178: INTEGER I, IERR, ITS, J
179: DOUBLE PRECISION GROWTO, NRMSML, ROOTN, RTEMP, SCALE, VNORM
180: COMPLEX*16 CDUM, EI, EJ, TEMP, X
181: * ..
182: * .. External Functions ..
183: INTEGER IZAMAX
184: DOUBLE PRECISION DZASUM, DZNRM2
185: COMPLEX*16 ZLADIV
186: EXTERNAL IZAMAX, DZASUM, DZNRM2, ZLADIV
187: * ..
188: * .. External Subroutines ..
189: EXTERNAL ZDSCAL, ZLATRS
190: * ..
191: * .. Intrinsic Functions ..
192: INTRINSIC ABS, DBLE, DIMAG, MAX, SQRT
193: * ..
194: * .. Statement Functions ..
195: DOUBLE PRECISION CABS1
196: * ..
197: * .. Statement Function definitions ..
198: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
199: * ..
200: * .. Executable Statements ..
201: *
202: INFO = 0
203: *
204: * GROWTO is the threshold used in the acceptance test for an
205: * eigenvector.
206: *
207: ROOTN = SQRT( DBLE( N ) )
208: GROWTO = TENTH / ROOTN
209: NRMSML = MAX( ONE, EPS3*ROOTN )*SMLNUM
210: *
211: * Form B = H - W*I (except that the subdiagonal elements are not
212: * stored).
213: *
214: DO 20 J = 1, N
215: DO 10 I = 1, J - 1
216: B( I, J ) = H( I, J )
217: 10 CONTINUE
218: B( J, J ) = H( J, J ) - W
219: 20 CONTINUE
220: *
221: IF( NOINIT ) THEN
222: *
223: * Initialize V.
224: *
225: DO 30 I = 1, N
226: V( I ) = EPS3
227: 30 CONTINUE
228: ELSE
229: *
230: * Scale supplied initial vector.
231: *
232: VNORM = DZNRM2( N, V, 1 )
233: CALL ZDSCAL( N, ( EPS3*ROOTN ) / MAX( VNORM, NRMSML ), V, 1 )
234: END IF
235: *
236: IF( RIGHTV ) THEN
237: *
238: * LU decomposition with partial pivoting of B, replacing zero
239: * pivots by EPS3.
240: *
241: DO 60 I = 1, N - 1
242: EI = H( I+1, I )
243: IF( CABS1( B( I, I ) ).LT.CABS1( EI ) ) THEN
244: *
245: * Interchange rows and eliminate.
246: *
247: X = ZLADIV( B( I, I ), EI )
248: B( I, I ) = EI
249: DO 40 J = I + 1, N
250: TEMP = B( I+1, J )
251: B( I+1, J ) = B( I, J ) - X*TEMP
252: B( I, J ) = TEMP
253: 40 CONTINUE
254: ELSE
255: *
256: * Eliminate without interchange.
257: *
258: IF( B( I, I ).EQ.ZERO )
259: $ B( I, I ) = EPS3
260: X = ZLADIV( EI, B( I, I ) )
261: IF( X.NE.ZERO ) THEN
262: DO 50 J = I + 1, N
263: B( I+1, J ) = B( I+1, J ) - X*B( I, J )
264: 50 CONTINUE
265: END IF
266: END IF
267: 60 CONTINUE
268: IF( B( N, N ).EQ.ZERO )
269: $ B( N, N ) = EPS3
270: *
271: TRANS = 'N'
272: *
273: ELSE
274: *
275: * UL decomposition with partial pivoting of B, replacing zero
276: * pivots by EPS3.
277: *
278: DO 90 J = N, 2, -1
279: EJ = H( J, J-1 )
280: IF( CABS1( B( J, J ) ).LT.CABS1( EJ ) ) THEN
281: *
282: * Interchange columns and eliminate.
283: *
284: X = ZLADIV( B( J, J ), EJ )
285: B( J, J ) = EJ
286: DO 70 I = 1, J - 1
287: TEMP = B( I, J-1 )
288: B( I, J-1 ) = B( I, J ) - X*TEMP
289: B( I, J ) = TEMP
290: 70 CONTINUE
291: ELSE
292: *
293: * Eliminate without interchange.
294: *
295: IF( B( J, J ).EQ.ZERO )
296: $ B( J, J ) = EPS3
297: X = ZLADIV( EJ, B( J, J ) )
298: IF( X.NE.ZERO ) THEN
299: DO 80 I = 1, J - 1
300: B( I, J-1 ) = B( I, J-1 ) - X*B( I, J )
301: 80 CONTINUE
302: END IF
303: END IF
304: 90 CONTINUE
305: IF( B( 1, 1 ).EQ.ZERO )
306: $ B( 1, 1 ) = EPS3
307: *
308: TRANS = 'C'
309: *
310: END IF
311: *
312: NORMIN = 'N'
313: DO 110 ITS = 1, N
314: *
315: * Solve U*x = scale*v for a right eigenvector
1.8 bertrand 316: * or U**H *x = scale*v for a left eigenvector,
1.1 bertrand 317: * overwriting x on v.
318: *
319: CALL ZLATRS( 'Upper', TRANS, 'Nonunit', NORMIN, N, B, LDB, V,
320: $ SCALE, RWORK, IERR )
321: NORMIN = 'Y'
322: *
323: * Test for sufficient growth in the norm of v.
324: *
325: VNORM = DZASUM( N, V, 1 )
326: IF( VNORM.GE.GROWTO*SCALE )
327: $ GO TO 120
328: *
329: * Choose new orthogonal starting vector and try again.
330: *
331: RTEMP = EPS3 / ( ROOTN+ONE )
332: V( 1 ) = EPS3
333: DO 100 I = 2, N
334: V( I ) = RTEMP
335: 100 CONTINUE
336: V( N-ITS+1 ) = V( N-ITS+1 ) - EPS3*ROOTN
337: 110 CONTINUE
338: *
339: * Failure to find eigenvector in N iterations.
340: *
341: INFO = 1
342: *
343: 120 CONTINUE
344: *
345: * Normalize eigenvector.
346: *
347: I = IZAMAX( N, V, 1 )
348: CALL ZDSCAL( N, ONE / CABS1( V( I ) ), V, 1 )
349: *
350: RETURN
351: *
352: * End of ZLAEIN
353: *
354: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlaein.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack