LAPACK  3.10.0
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
9 *> Download ZHEEV + dependencies
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 *> \ingroup complex16HEeigen
136 *
137 * =====================================================================
138  SUBROUTINE zheev( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
139  $ INFO )
140 *
141 * -- LAPACK driver routine --
142 * -- LAPACK is a software package provided by Univ. of Tennessee, --
143 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
144 *
145 * .. Scalar Arguments ..
146  CHARACTER JOBZ, UPLO
147  INTEGER INFO, LDA, LWORK, N
148 * ..
149 * .. Array Arguments ..
150  DOUBLE PRECISION RWORK( * ), W( * )
151  COMPLEX*16 A( LDA, * ), WORK( * )
152 * ..
153 *
154 * =====================================================================
155 *
156 * .. Parameters ..
157  DOUBLE PRECISION ZERO, ONE
158  parameter( zero = 0.0d0, one = 1.0d0 )
159  COMPLEX*16 CONE
160  parameter( cone = ( 1.0d0, 0.0d0 ) )
161 * ..
162 * .. Local Scalars ..
163  LOGICAL LOWER, LQUERY, WANTZ
164  INTEGER IINFO, IMAX, INDE, INDTAU, INDWRK, ISCALE,
165  $ llwork, lwkopt, nb
166  DOUBLE PRECISION ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA,
167  $ smlnum
168 * ..
169 * .. External Functions ..
170  LOGICAL LSAME
171  INTEGER ILAENV
172  DOUBLE PRECISION DLAMCH, ZLANHE
173  EXTERNAL lsame, ilaenv, dlamch, zlanhe
174 * ..
175 * .. External Subroutines ..
176  EXTERNAL dscal, dsterf, xerbla, zhetrd, zlascl, zsteqr,
177  $ zungtr
178 * ..
179 * .. Intrinsic Functions ..
180  INTRINSIC max, sqrt
181 * ..
182 * .. Executable Statements ..
183 *
184 * Test the input parameters.
185 *
186  wantz = lsame( jobz, 'V' )
187  lower = lsame( uplo, 'L' )
188  lquery = ( lwork.EQ.-1 )
189 *
190  info = 0
191  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
192  info = -1
193  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
194  info = -2
195  ELSE IF( n.LT.0 ) THEN
196  info = -3
197  ELSE IF( lda.LT.max( 1, n ) ) THEN
198  info = -5
199  END IF
200 *
201  IF( info.EQ.0 ) THEN
202  nb = ilaenv( 1, 'ZHETRD', uplo, n, -1, -1, -1 )
203  lwkopt = max( 1, ( nb+1 )*n )
204  work( 1 ) = lwkopt
205 *
206  IF( lwork.LT.max( 1, 2*n-1 ) .AND. .NOT.lquery )
207  $ info = -8
208  END IF
209 *
210  IF( info.NE.0 ) THEN
211  CALL xerbla( 'ZHEEV ', -info )
212  RETURN
213  ELSE IF( lquery ) THEN
214  RETURN
215  END IF
216 *
217 * Quick return if possible
218 *
219  IF( n.EQ.0 ) THEN
220  RETURN
221  END IF
222 *
223  IF( n.EQ.1 ) THEN
224  w( 1 ) = dble( a( 1, 1 ) )
225  work( 1 ) = 1
226  IF( wantz )
227  $ a( 1, 1 ) = cone
228  RETURN
229  END IF
230 *
231 * Get machine constants.
232 *
233  safmin = dlamch( 'Safe minimum' )
234  eps = dlamch( 'Precision' )
235  smlnum = safmin / eps
236  bignum = one / smlnum
237  rmin = sqrt( smlnum )
238  rmax = sqrt( bignum )
239 *
240 * Scale matrix to allowable range, if necessary.
241 *
242  anrm = zlanhe( 'M', uplo, n, a, lda, rwork )
243  iscale = 0
244  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
245  iscale = 1
246  sigma = rmin / anrm
247  ELSE IF( anrm.GT.rmax ) THEN
248  iscale = 1
249  sigma = rmax / anrm
250  END IF
251  IF( iscale.EQ.1 )
252  $ CALL zlascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
253 *
254 * Call ZHETRD to reduce Hermitian matrix to tridiagonal form.
255 *
256  inde = 1
257  indtau = 1
258  indwrk = indtau + n
259  llwork = lwork - indwrk + 1
260  CALL zhetrd( uplo, n, a, lda, w, rwork( inde ), work( indtau ),
261  $ work( indwrk ), llwork, iinfo )
262 *
263 * For eigenvalues only, call DSTERF. For eigenvectors, first call
264 * ZUNGTR to generate the unitary matrix, then call ZSTEQR.
265 *
266  IF( .NOT.wantz ) THEN
267  CALL dsterf( n, w, rwork( inde ), info )
268  ELSE
269  CALL zungtr( uplo, n, a, lda, work( indtau ), work( indwrk ),
270  $ llwork, iinfo )
271  indwrk = inde + n
272  CALL zsteqr( jobz, n, w, rwork( inde ), a, lda,
273  $ rwork( indwrk ), info )
274  END IF
275 *
276 * If matrix was scaled, then rescale eigenvalues appropriately.
277 *
278  IF( iscale.EQ.1 ) THEN
279  IF( info.EQ.0 ) THEN
280  imax = n
281  ELSE
282  imax = info - 1
283  END IF
284  CALL dscal( imax, one / sigma, w, 1 )
285  END IF
286 *
287 * Set WORK(1) to optimal complex workspace size.
288 *
289  work( 1 ) = lwkopt
290 *
291  RETURN
292 *
293 * End of ZHEEV
294 *
295  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:86
subroutine zhetrd(UPLO, N, A, LDA, D, E, TAU, WORK, LWORK, INFO)
ZHETRD
Definition: zhetrd.f:192
subroutine zheev(JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO)
ZHEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices
Definition: zheev.f:140
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:143
subroutine zsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
ZSTEQR
Definition: zsteqr.f:132
subroutine zungtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGTR
Definition: zungtr.f:123
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79