LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
zheev.f
Go to the documentation of this file.
1 *> \brief <b> ZHEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices</b>
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zheev.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zheev.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zheev.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
22 * INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER JOBZ, UPLO
26 * INTEGER INFO, LDA, LWORK, N
27 * ..
28 * .. Array Arguments ..
29 * DOUBLE PRECISION RWORK( * ), W( * )
30 * COMPLEX*16 A( LDA, * ), WORK( * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> ZHEEV computes all eigenvalues and, optionally, eigenvectors of a
40 *> complex Hermitian matrix A.
41 *> \endverbatim
42 *
43 * Arguments:
44 * ==========
45 *
46 *> \param[in] JOBZ
47 *> \verbatim
48 *> JOBZ is CHARACTER*1
49 *> = 'N': Compute eigenvalues only;
50 *> = 'V': Compute eigenvalues and eigenvectors.
51 *> \endverbatim
52 *>
53 *> \param[in] UPLO
54 *> \verbatim
55 *> UPLO is CHARACTER*1
56 *> = 'U': Upper triangle of A is stored;
57 *> = 'L': Lower triangle of A is stored.
58 *> \endverbatim
59 *>
60 *> \param[in] N
61 *> \verbatim
62 *> N is INTEGER
63 *> The order of the matrix A. N >= 0.
64 *> \endverbatim
65 *>
66 *> \param[in,out] A
67 *> \verbatim
68 *> A is COMPLEX*16 array, dimension (LDA, N)
69 *> On entry, the Hermitian matrix A. If UPLO = 'U', the
70 *> leading N-by-N upper triangular part of A contains the
71 *> upper triangular part of the matrix A. If UPLO = 'L',
72 *> the leading N-by-N lower triangular part of A contains
73 *> the lower triangular part of the matrix A.
74 *> On exit, if JOBZ = 'V', then if INFO = 0, A contains the
75 *> orthonormal eigenvectors of the matrix A.
76 *> If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
77 *> or the upper triangle (if UPLO='U') of A, including the
78 *> diagonal, is destroyed.
79 *> \endverbatim
80 *>
81 *> \param[in] LDA
82 *> \verbatim
83 *> LDA is INTEGER
84 *> The leading dimension of the array A. LDA >= max(1,N).
85 *> \endverbatim
86 *>
87 *> \param[out] W
88 *> \verbatim
89 *> W is DOUBLE PRECISION array, dimension (N)
90 *> If INFO = 0, the eigenvalues in ascending order.
91 *> \endverbatim
92 *>
93 *> \param[out] WORK
94 *> \verbatim
95 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
96 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
97 *> \endverbatim
98 *>
99 *> \param[in] LWORK
100 *> \verbatim
101 *> LWORK is INTEGER
102 *> The length of the array WORK. LWORK >= max(1,2*N-1).
103 *> For optimal efficiency, LWORK >= (NB+1)*N,
104 *> where NB is the blocksize for ZHETRD returned by ILAENV.
105 *>
106 *> If LWORK = -1, then a workspace query is assumed; the routine
107 *> only calculates the optimal size of the WORK array, returns
108 *> this value as the first entry of the WORK array, and no error
109 *> message related to LWORK is issued by XERBLA.
110 *> \endverbatim
111 *>
112 *> \param[out] RWORK
113 *> \verbatim
114 *> RWORK is DOUBLE PRECISION array, dimension (max(1, 3*N-2))
115 *> \endverbatim
116 *>
117 *> \param[out] INFO
118 *> \verbatim
119 *> INFO is INTEGER
120 *> = 0: successful exit
121 *> < 0: if INFO = -i, the i-th argument had an illegal value
122 *> > 0: if INFO = i, the algorithm failed to converge; i
123 *> off-diagonal elements of an intermediate tridiagonal
124 *> form did not converge to zero.
125 *> \endverbatim
126 *
127 * Authors:
128 * ========
129 *
130 *> \author Univ. of Tennessee
131 *> \author Univ. of California Berkeley
132 *> \author Univ. of Colorado Denver
133 *> \author NAG Ltd.
134 *
135 *> \date November 2011
136 *
137 *> \ingroup complex16HEeigen
138 *
139 * =====================================================================
140  SUBROUTINE zheev( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
141  \$ info )
142 *
143 * -- LAPACK driver routine (version 3.4.0) --
144 * -- LAPACK is a software package provided by Univ. of Tennessee, --
145 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
146 * November 2011
147 *
148 * .. Scalar Arguments ..
149  CHARACTER jobz, uplo
150  INTEGER info, lda, lwork, n
151 * ..
152 * .. Array Arguments ..
153  DOUBLE PRECISION rwork( * ), w( * )
154  COMPLEX*16 a( lda, * ), work( * )
155 * ..
156 *
157 * =====================================================================
158 *
159 * .. Parameters ..
160  DOUBLE PRECISION zero, one
161  parameter( zero = 0.0d0, one = 1.0d0 )
162  COMPLEX*16 cone
163  parameter( cone = ( 1.0d0, 0.0d0 ) )
164 * ..
165 * .. Local Scalars ..
166  LOGICAL lower, lquery, wantz
167  INTEGER iinfo, imax, inde, indtau, indwrk, iscale,
168  \$ llwork, lwkopt, nb
169  DOUBLE PRECISION anrm, bignum, eps, rmax, rmin, safmin, sigma,
170  \$ smlnum
171 * ..
172 * .. External Functions ..
173  LOGICAL lsame
174  INTEGER ilaenv
175  DOUBLE PRECISION dlamch, zlanhe
176  EXTERNAL lsame, ilaenv, dlamch, zlanhe
177 * ..
178 * .. External Subroutines ..
179  EXTERNAL dscal, dsterf, xerbla, zhetrd, zlascl, zsteqr,
180  \$ zungtr
181 * ..
182 * .. Intrinsic Functions ..
183  INTRINSIC max, sqrt
184 * ..
185 * .. Executable Statements ..
186 *
187 * Test the input parameters.
188 *
189  wantz = lsame( jobz, 'V' )
190  lower = lsame( uplo, 'L' )
191  lquery = ( lwork.EQ.-1 )
192 *
193  info = 0
194  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
195  info = -1
196  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
197  info = -2
198  ELSE IF( n.LT.0 ) THEN
199  info = -3
200  ELSE IF( lda.LT.max( 1, n ) ) THEN
201  info = -5
202  END IF
203 *
204  IF( info.EQ.0 ) THEN
205  nb = ilaenv( 1, 'ZHETRD', uplo, n, -1, -1, -1 )
206  lwkopt = max( 1, ( nb+1 )*n )
207  work( 1 ) = lwkopt
208 *
209  IF( lwork.LT.max( 1, 2*n-1 ) .AND. .NOT.lquery )
210  \$ info = -8
211  END IF
212 *
213  IF( info.NE.0 ) THEN
214  CALL xerbla( 'ZHEEV ', -info )
215  return
216  ELSE IF( lquery ) THEN
217  return
218  END IF
219 *
220 * Quick return if possible
221 *
222  IF( n.EQ.0 ) THEN
223  return
224  END IF
225 *
226  IF( n.EQ.1 ) THEN
227  w( 1 ) = a( 1, 1 )
228  work( 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  llwork = lwork - indwrk + 1
263  CALL zhetrd( uplo, n, a, lda, w, rwork( inde ), work( indtau ),
264  \$ work( indwrk ), llwork, iinfo )
265 *
266 * For eigenvalues only, call DSTERF. For eigenvectors, first call
267 * ZUNGTR to generate the unitary matrix, then call ZSTEQR.
268 *
269  IF( .NOT.wantz ) THEN
270  CALL dsterf( n, w, rwork( inde ), info )
271  ELSE
272  CALL zungtr( uplo, n, a, lda, work( indtau ), work( indwrk ),
273  \$ llwork, iinfo )
274  indwrk = inde + n
275  CALL zsteqr( jobz, n, w, rwork( inde ), a, lda,
276  \$ rwork( indwrk ), info )
277  END IF
278 *
279 * If matrix was scaled, then rescale eigenvalues appropriately.
280 *
281  IF( iscale.EQ.1 ) THEN
282  IF( info.EQ.0 ) THEN
283  imax = n
284  ELSE
285  imax = info - 1
286  END IF
287  CALL dscal( imax, one / sigma, w, 1 )
288  END IF
289 *
290 * Set WORK(1) to optimal complex workspace size.
291 *
292  work( 1 ) = lwkopt
293 *
294  return
295 *
296 * End of ZHEEV
297 *
298  END