1: SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
2: $ INFO )
3: *
4: * -- LAPACK driver routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: CHARACTER JOBZ, UPLO
11: INTEGER INFO, LDA, LWORK, N
12: * ..
13: * .. Array Arguments ..
14: DOUBLE PRECISION RWORK( * ), W( * )
15: COMPLEX*16 A( LDA, * ), WORK( * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * ZHEEV computes all eigenvalues and, optionally, eigenvectors of a
22: * complex Hermitian matrix A.
23: *
24: * Arguments
25: * =========
26: *
27: * JOBZ (input) CHARACTER*1
28: * = 'N': Compute eigenvalues only;
29: * = 'V': Compute eigenvalues and eigenvectors.
30: *
31: * UPLO (input) CHARACTER*1
32: * = 'U': Upper triangle of A is stored;
33: * = 'L': Lower triangle of A is stored.
34: *
35: * N (input) INTEGER
36: * The order of the matrix A. N >= 0.
37: *
38: * A (input/output) COMPLEX*16 array, dimension (LDA, N)
39: * On entry, the Hermitian matrix A. If UPLO = 'U', the
40: * leading N-by-N upper triangular part of A contains the
41: * upper triangular part of the matrix A. If UPLO = 'L',
42: * the leading N-by-N lower triangular part of A contains
43: * the lower triangular part of the matrix A.
44: * On exit, if JOBZ = 'V', then if INFO = 0, A contains the
45: * orthonormal eigenvectors of the matrix A.
46: * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
47: * or the upper triangle (if UPLO='U') of A, including the
48: * diagonal, is destroyed.
49: *
50: * LDA (input) INTEGER
51: * The leading dimension of the array A. LDA >= max(1,N).
52: *
53: * W (output) DOUBLE PRECISION array, dimension (N)
54: * If INFO = 0, the eigenvalues in ascending order.
55: *
56: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
57: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
58: *
59: * LWORK (input) INTEGER
60: * The length of the array WORK. LWORK >= max(1,2*N-1).
61: * For optimal efficiency, LWORK >= (NB+1)*N,
62: * where NB is the blocksize for ZHETRD returned by ILAENV.
63: *
64: * If LWORK = -1, then a workspace query is assumed; the routine
65: * only calculates the optimal size of the WORK array, returns
66: * this value as the first entry of the WORK array, and no error
67: * message related to LWORK is issued by XERBLA.
68: *
69: * RWORK (workspace) DOUBLE PRECISION array, dimension (max(1, 3*N-2))
70: *
71: * INFO (output) INTEGER
72: * = 0: successful exit
73: * < 0: if INFO = -i, the i-th argument had an illegal value
74: * > 0: if INFO = i, the algorithm failed to converge; i
75: * off-diagonal elements of an intermediate tridiagonal
76: * form did not converge to zero.
77: *
78: * =====================================================================
79: *
80: * .. Parameters ..
81: DOUBLE PRECISION ZERO, ONE
82: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
83: COMPLEX*16 CONE
84: PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) )
85: * ..
86: * .. Local Scalars ..
87: LOGICAL LOWER, LQUERY, WANTZ
88: INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
89: $ LLWORK, LWKOPT, NB
90: DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
91: $ SMLNUM
92: * ..
93: * .. External Functions ..
94: LOGICAL LSAME
95: INTEGER ILAENV
96: DOUBLE PRECISION DLAMCH, ZLANHE
97: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANHE
98: * ..
99: * .. External Subroutines ..
100: EXTERNAL DSCAL, DSTERF, XERBLA, ZHETRD, ZLASCL, ZSTEQR,
101: $ ZUNGTR
102: * ..
103: * .. Intrinsic Functions ..
104: INTRINSIC MAX, SQRT
105: * ..
106: * .. Executable Statements ..
107: *
108: * Test the input parameters.
109: *
110: WANTZ = LSAME( JOBZ, 'V' )
111: LOWER = LSAME( UPLO, 'L' )
112: LQUERY = ( LWORK.EQ.-1 )
113: *
114: INFO = 0
115: IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
116: INFO = -1
117: ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
118: INFO = -2
119: ELSE IF( N.LT.0 ) THEN
120: INFO = -3
121: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
122: INFO = -5
123: END IF
124: *
125: IF( INFO.EQ.0 ) THEN
126: NB = ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 )
127: LWKOPT = MAX( 1, ( NB+1 )*N )
128: WORK( 1 ) = LWKOPT
129: *
130: IF( LWORK.LT.MAX( 1, 2*N-1 ) .AND. .NOT.LQUERY )
131: $ INFO = -8
132: END IF
133: *
134: IF( INFO.NE.0 ) THEN
135: CALL XERBLA( 'ZHEEV ', -INFO )
136: RETURN
137: ELSE IF( LQUERY ) THEN
138: RETURN
139: END IF
140: *
141: * Quick return if possible
142: *
143: IF( N.EQ.0 ) THEN
144: RETURN
145: END IF
146: *
147: IF( N.EQ.1 ) THEN
148: W( 1 ) = A( 1, 1 )
149: WORK( 1 ) = 1
150: IF( WANTZ )
151: $ A( 1, 1 ) = CONE
152: RETURN
153: END IF
154: *
155: * Get machine constants.
156: *
157: SAFMIN = DLAMCH( 'Safe minimum' )
158: EPS = DLAMCH( 'Precision' )
159: SMLNUM = SAFMIN / EPS
160: BIGNUM = ONE / SMLNUM
161: RMIN = SQRT( SMLNUM )
162: RMAX = SQRT( BIGNUM )
163: *
164: * Scale matrix to allowable range, if necessary.
165: *
166: ANRM = ZLANHE( 'M', UPLO, N, A, LDA, RWORK )
167: ISCALE = 0
168: IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
169: ISCALE = 1
170: SIGMA = RMIN / ANRM
171: ELSE IF( ANRM.GT.RMAX ) THEN
172: ISCALE = 1
173: SIGMA = RMAX / ANRM
174: END IF
175: IF( ISCALE.EQ.1 )
176: $ CALL ZLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
177: *
178: * Call ZHETRD to reduce Hermitian matrix to tridiagonal form.
179: *
180: INDE = 1
181: INDTAU = 1
182: INDWRK = INDTAU + N
183: LLWORK = LWORK - INDWRK + 1
184: CALL ZHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ),
185: $ WORK( INDWRK ), LLWORK, IINFO )
186: *
187: * For eigenvalues only, call DSTERF. For eigenvectors, first call
188: * ZUNGTR to generate the unitary matrix, then call ZSTEQR.
189: *
190: IF( .NOT.WANTZ ) THEN
191: CALL DSTERF( N, W, RWORK( INDE ), INFO )
192: ELSE
193: CALL ZUNGTR( UPLO, N, A, LDA, WORK( INDTAU ), WORK( INDWRK ),
194: $ LLWORK, IINFO )
195: INDWRK = INDE + N
196: CALL ZSTEQR( JOBZ, N, W, RWORK( INDE ), A, LDA,
197: $ RWORK( INDWRK ), INFO )
198: END IF
199: *
200: * If matrix was scaled, then rescale eigenvalues appropriately.
201: *
202: IF( ISCALE.EQ.1 ) THEN
203: IF( INFO.EQ.0 ) THEN
204: IMAX = N
205: ELSE
206: IMAX = INFO - 1
207: END IF
208: CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
209: END IF
210: *
211: * Set WORK(1) to optimal complex workspace size.
212: *
213: WORK( 1 ) = LWKOPT
214: *
215: RETURN
216: *
217: * End of ZHEEV
218: *
219: END
CVSweb interface <joel.bertrand@systella.fr>