Annotation of rpl/lapack/lapack/zheevd.f, revision 1.4
1.1 bertrand 1: SUBROUTINE ZHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
2: $ LRWORK, IWORK, LIWORK, 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, LIWORK, LRWORK, LWORK, N
12: * ..
13: * .. Array Arguments ..
14: INTEGER IWORK( * )
15: DOUBLE PRECISION RWORK( * ), W( * )
16: COMPLEX*16 A( LDA, * ), WORK( * )
17: * ..
18: *
19: * Purpose
20: * =======
21: *
22: * ZHEEVD computes all eigenvalues and, optionally, eigenvectors of a
23: * complex Hermitian matrix A. If eigenvectors are desired, it uses a
24: * divide and conquer algorithm.
25: *
26: * The divide and conquer algorithm makes very mild assumptions about
27: * floating point arithmetic. It will work on machines with a guard
28: * digit in add/subtract, or on those binary machines without guard
29: * digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
30: * Cray-2. It could conceivably fail on hexadecimal or decimal machines
31: * without guard digits, but we know of none.
32: *
33: * Arguments
34: * =========
35: *
36: * JOBZ (input) CHARACTER*1
37: * = 'N': Compute eigenvalues only;
38: * = 'V': Compute eigenvalues and eigenvectors.
39: *
40: * UPLO (input) CHARACTER*1
41: * = 'U': Upper triangle of A is stored;
42: * = 'L': Lower triangle of A is stored.
43: *
44: * N (input) INTEGER
45: * The order of the matrix A. N >= 0.
46: *
47: * A (input/output) COMPLEX*16 array, dimension (LDA, N)
48: * On entry, the Hermitian matrix A. If UPLO = 'U', the
49: * leading N-by-N upper triangular part of A contains the
50: * upper triangular part of the matrix A. If UPLO = 'L',
51: * the leading N-by-N lower triangular part of A contains
52: * the lower triangular part of the matrix A.
53: * On exit, if JOBZ = 'V', then if INFO = 0, A contains the
54: * orthonormal eigenvectors of the matrix A.
55: * If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
56: * or the upper triangle (if UPLO='U') of A, including the
57: * diagonal, is destroyed.
58: *
59: * LDA (input) INTEGER
60: * The leading dimension of the array A. LDA >= max(1,N).
61: *
62: * W (output) DOUBLE PRECISION array, dimension (N)
63: * If INFO = 0, the eigenvalues in ascending order.
64: *
65: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
66: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
67: *
68: * LWORK (input) INTEGER
69: * The length of the array WORK.
70: * If N <= 1, LWORK must be at least 1.
71: * If JOBZ = 'N' and N > 1, LWORK must be at least N + 1.
72: * If JOBZ = 'V' and N > 1, LWORK must be at least 2*N + N**2.
73: *
74: * If LWORK = -1, then a workspace query is assumed; the routine
75: * only calculates the optimal sizes of the WORK, RWORK and
76: * IWORK arrays, returns these values as the first entries of
77: * the WORK, RWORK and IWORK arrays, and no error message
78: * related to LWORK or LRWORK or LIWORK is issued by XERBLA.
79: *
80: * RWORK (workspace/output) DOUBLE PRECISION array,
81: * dimension (LRWORK)
82: * On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
83: *
84: * LRWORK (input) INTEGER
85: * The dimension of the array RWORK.
86: * If N <= 1, LRWORK must be at least 1.
87: * If JOBZ = 'N' and N > 1, LRWORK must be at least N.
88: * If JOBZ = 'V' and N > 1, LRWORK must be at least
89: * 1 + 5*N + 2*N**2.
90: *
91: * If LRWORK = -1, then a workspace query is assumed; the
92: * routine only calculates the optimal sizes of the WORK, RWORK
93: * and IWORK arrays, returns these values as the first entries
94: * of the WORK, RWORK and IWORK arrays, and no error message
95: * related to LWORK or LRWORK or LIWORK is issued by XERBLA.
96: *
97: * IWORK (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
98: * On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
99: *
100: * LIWORK (input) INTEGER
101: * The dimension of the array IWORK.
102: * If N <= 1, LIWORK must be at least 1.
103: * If JOBZ = 'N' and N > 1, LIWORK must be at least 1.
104: * If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
105: *
106: * If LIWORK = -1, then a workspace query is assumed; the
107: * routine only calculates the optimal sizes of the WORK, RWORK
108: * and IWORK arrays, returns these values as the first entries
109: * of the WORK, RWORK and IWORK arrays, and no error message
110: * related to LWORK or LRWORK or LIWORK is issued by XERBLA.
111: *
112: * INFO (output) INTEGER
113: * = 0: successful exit
114: * < 0: if INFO = -i, the i-th argument had an illegal value
115: * > 0: if INFO = i and JOBZ = 'N', then the algorithm failed
116: * to converge; i off-diagonal elements of an intermediate
117: * tridiagonal form did not converge to zero;
118: * if INFO = i and JOBZ = 'V', then the algorithm failed
119: * to compute an eigenvalue while working on the submatrix
120: * lying in rows and columns INFO/(N+1) through
121: * mod(INFO,N+1).
122: *
123: * Further Details
124: * ===============
125: *
126: * Based on contributions by
127: * Jeff Rutter, Computer Science Division, University of California
128: * at Berkeley, USA
129: *
130: * Modified description of INFO. Sven, 16 Feb 05.
131: * =====================================================================
132: *
133: * .. Parameters ..
134: DOUBLE PRECISION ZERO, ONE
135: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
136: COMPLEX*16 CONE
137: PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) )
138: * ..
139: * .. Local Scalars ..
140: LOGICAL LOWER, LQUERY, WANTZ
141: INTEGER IINFO, IMAX, INDE, INDRWK, INDTAU, INDWK2,
142: $ INDWRK, ISCALE, LIOPT, LIWMIN, LLRWK, LLWORK,
143: $ LLWRK2, LOPT, LROPT, LRWMIN, LWMIN
144: DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
145: $ SMLNUM
146: * ..
147: * .. External Functions ..
148: LOGICAL LSAME
149: INTEGER ILAENV
150: DOUBLE PRECISION DLAMCH, ZLANHE
151: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANHE
152: * ..
153: * .. External Subroutines ..
154: EXTERNAL DSCAL, DSTERF, XERBLA, ZHETRD, ZLACPY, ZLASCL,
155: $ ZSTEDC, ZUNMTR
156: * ..
157: * .. Intrinsic Functions ..
158: INTRINSIC MAX, SQRT
159: * ..
160: * .. Executable Statements ..
161: *
162: * Test the input parameters.
163: *
164: WANTZ = LSAME( JOBZ, 'V' )
165: LOWER = LSAME( UPLO, 'L' )
166: LQUERY = ( LWORK.EQ.-1 .OR. LRWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
167: *
168: INFO = 0
169: IF( .NOT.( WANTZ .OR. LSAME( JOBZ, 'N' ) ) ) THEN
170: INFO = -1
171: ELSE IF( .NOT.( LOWER .OR. LSAME( UPLO, 'U' ) ) ) THEN
172: INFO = -2
173: ELSE IF( N.LT.0 ) THEN
174: INFO = -3
175: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
176: INFO = -5
177: END IF
178: *
179: IF( INFO.EQ.0 ) THEN
180: IF( N.LE.1 ) THEN
181: LWMIN = 1
182: LRWMIN = 1
183: LIWMIN = 1
184: LOPT = LWMIN
185: LROPT = LRWMIN
186: LIOPT = LIWMIN
187: ELSE
188: IF( WANTZ ) THEN
189: LWMIN = 2*N + N*N
190: LRWMIN = 1 + 5*N + 2*N**2
191: LIWMIN = 3 + 5*N
192: ELSE
193: LWMIN = N + 1
194: LRWMIN = N
195: LIWMIN = 1
196: END IF
197: LOPT = MAX( LWMIN, N +
198: $ ILAENV( 1, 'ZHETRD', UPLO, N, -1, -1, -1 ) )
199: LROPT = LRWMIN
200: LIOPT = LIWMIN
201: END IF
202: WORK( 1 ) = LOPT
203: RWORK( 1 ) = LROPT
204: IWORK( 1 ) = LIOPT
205: *
206: IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
207: INFO = -8
208: ELSE IF( LRWORK.LT.LRWMIN .AND. .NOT.LQUERY ) THEN
209: INFO = -10
210: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
211: INFO = -12
212: END IF
213: END IF
214: *
215: IF( INFO.NE.0 ) THEN
216: CALL XERBLA( 'ZHEEVD', -INFO )
217: RETURN
218: ELSE IF( LQUERY ) THEN
219: RETURN
220: END IF
221: *
222: * Quick return if possible
223: *
224: IF( N.EQ.0 )
225: $ RETURN
226: *
227: IF( N.EQ.1 ) THEN
228: W( 1 ) = A( 1, 1 )
229: IF( WANTZ )
230: $ A( 1, 1 ) = CONE
231: RETURN
232: END IF
233: *
234: * Get machine constants.
235: *
236: SAFMIN = DLAMCH( 'Safe minimum' )
237: EPS = DLAMCH( 'Precision' )
238: SMLNUM = SAFMIN / EPS
239: BIGNUM = ONE / SMLNUM
240: RMIN = SQRT( SMLNUM )
241: RMAX = SQRT( BIGNUM )
242: *
243: * Scale matrix to allowable range, if necessary.
244: *
245: ANRM = ZLANHE( 'M', UPLO, N, A, LDA, RWORK )
246: ISCALE = 0
247: IF( ANRM.GT.ZERO .AND. ANRM.LT.RMIN ) THEN
248: ISCALE = 1
249: SIGMA = RMIN / ANRM
250: ELSE IF( ANRM.GT.RMAX ) THEN
251: ISCALE = 1
252: SIGMA = RMAX / ANRM
253: END IF
254: IF( ISCALE.EQ.1 )
255: $ CALL ZLASCL( UPLO, 0, 0, ONE, SIGMA, N, N, A, LDA, INFO )
256: *
257: * Call ZHETRD to reduce Hermitian matrix to tridiagonal form.
258: *
259: INDE = 1
260: INDTAU = 1
261: INDWRK = INDTAU + N
262: INDRWK = INDE + N
263: INDWK2 = INDWRK + N*N
264: LLWORK = LWORK - INDWRK + 1
265: LLWRK2 = LWORK - INDWK2 + 1
266: LLRWK = LRWORK - INDRWK + 1
267: CALL ZHETRD( UPLO, N, A, LDA, W, RWORK( INDE ), WORK( INDTAU ),
268: $ WORK( INDWRK ), LLWORK, IINFO )
269: *
270: * For eigenvalues only, call DSTERF. For eigenvectors, first call
271: * ZSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
272: * tridiagonal matrix, then call ZUNMTR to multiply it to the
273: * Householder transformations represented as Householder vectors in
274: * A.
275: *
276: IF( .NOT.WANTZ ) THEN
277: CALL DSTERF( N, W, RWORK( INDE ), INFO )
278: ELSE
279: CALL ZSTEDC( 'I', N, W, RWORK( INDE ), WORK( INDWRK ), N,
280: $ WORK( INDWK2 ), LLWRK2, RWORK( INDRWK ), LLRWK,
281: $ IWORK, LIWORK, INFO )
282: CALL ZUNMTR( 'L', UPLO, 'N', N, N, A, LDA, WORK( INDTAU ),
283: $ WORK( INDWRK ), N, WORK( INDWK2 ), LLWRK2, IINFO )
284: CALL ZLACPY( 'A', N, N, WORK( INDWRK ), N, A, LDA )
285: END IF
286: *
287: * If matrix was scaled, then rescale eigenvalues appropriately.
288: *
289: IF( ISCALE.EQ.1 ) THEN
290: IF( INFO.EQ.0 ) THEN
291: IMAX = N
292: ELSE
293: IMAX = INFO - 1
294: END IF
295: CALL DSCAL( IMAX, ONE / SIGMA, W, 1 )
296: END IF
297: *
298: WORK( 1 ) = LOPT
299: RWORK( 1 ) = LROPT
300: IWORK( 1 ) = LIOPT
301: *
302: RETURN
303: *
304: * End of ZHEEVD
305: *
306: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zheevd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack