LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zheevd.f
Go to the documentation of this file.
1 *> \brief <b> ZHEEVD 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 ZHEEVD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zheevd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zheevd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zheevd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHEEVD( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
22 * LRWORK, IWORK, LIWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER JOBZ, UPLO
26 * INTEGER INFO, LDA, LIWORK, LRWORK, LWORK, N
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IWORK( * )
30 * DOUBLE PRECISION RWORK( * ), W( * )
31 * COMPLEX*16 A( LDA, * ), WORK( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> ZHEEVD computes all eigenvalues and, optionally, eigenvectors of a
41 *> complex Hermitian matrix A. If eigenvectors are desired, it uses a
42 *> divide and conquer algorithm.
43 *>
44 *> The divide and conquer algorithm makes very mild assumptions about
45 *> floating point arithmetic. It will work on machines with a guard
46 *> digit in add/subtract, or on those binary machines without guard
47 *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
48 *> Cray-2. It could conceivably fail on hexadecimal or decimal machines
49 *> without guard digits, but we know of none.
50 *> \endverbatim
51 *
52 * Arguments:
53 * ==========
54 *
55 *> \param[in] JOBZ
56 *> \verbatim
57 *> JOBZ is CHARACTER*1
58 *> = 'N': Compute eigenvalues only;
59 *> = 'V': Compute eigenvalues and eigenvectors.
60 *> \endverbatim
61 *>
62 *> \param[in] UPLO
63 *> \verbatim
64 *> UPLO is CHARACTER*1
65 *> = 'U': Upper triangle of A is stored;
66 *> = 'L': Lower triangle of A is stored.
67 *> \endverbatim
68 *>
69 *> \param[in] N
70 *> \verbatim
71 *> N is INTEGER
72 *> The order of the matrix A. N >= 0.
73 *> \endverbatim
74 *>
75 *> \param[in,out] A
76 *> \verbatim
77 *> A is COMPLEX*16 array, dimension (LDA, N)
78 *> On entry, the Hermitian matrix A. If UPLO = 'U', the
79 *> leading N-by-N upper triangular part of A contains the
80 *> upper triangular part of the matrix A. If UPLO = 'L',
81 *> the leading N-by-N lower triangular part of A contains
82 *> the lower triangular part of the matrix A.
83 *> On exit, if JOBZ = 'V', then if INFO = 0, A contains the
84 *> orthonormal eigenvectors of the matrix A.
85 *> If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
86 *> or the upper triangle (if UPLO='U') of A, including the
87 *> diagonal, is destroyed.
88 *> \endverbatim
89 *>
90 *> \param[in] LDA
91 *> \verbatim
92 *> LDA is INTEGER
93 *> The leading dimension of the array A. LDA >= max(1,N).
94 *> \endverbatim
95 *>
96 *> \param[out] W
97 *> \verbatim
98 *> W is DOUBLE PRECISION array, dimension (N)
99 *> If INFO = 0, the eigenvalues in ascending order.
100 *> \endverbatim
101 *>
102 *> \param[out] WORK
103 *> \verbatim
104 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
105 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
106 *> \endverbatim
107 *>
108 *> \param[in] LWORK
109 *> \verbatim
110 *> LWORK is INTEGER
111 *> The length of the array WORK.
112 *> If N <= 1, LWORK must be at least 1.
113 *> If JOBZ = 'N' and N > 1, LWORK must be at least N + 1.
114 *> If JOBZ = 'V' and N > 1, LWORK must be at least 2*N + N**2.
115 *>
116 *> If LWORK = -1, then a workspace query is assumed; the routine
117 *> only calculates the optimal sizes of the WORK, RWORK and
118 *> IWORK arrays, returns these values as the first entries of
119 *> the WORK, RWORK and IWORK arrays, and no error message
120 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
121 *> \endverbatim
122 *>
123 *> \param[out] RWORK
124 *> \verbatim
125 *> RWORK is DOUBLE PRECISION array,
126 *> dimension (LRWORK)
127 *> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
128 *> \endverbatim
129 *>
130 *> \param[in] LRWORK
131 *> \verbatim
132 *> LRWORK is INTEGER
133 *> The dimension of the array RWORK.
134 *> If N <= 1, LRWORK must be at least 1.
135 *> If JOBZ = 'N' and N > 1, LRWORK must be at least N.
136 *> If JOBZ = 'V' and N > 1, LRWORK must be at least
137 *> 1 + 5*N + 2*N**2.
138 *>
139 *> If LRWORK = -1, then a workspace query is assumed; the
140 *> routine only calculates the optimal sizes of the WORK, RWORK
141 *> and IWORK arrays, returns these values as the first entries
142 *> of the WORK, RWORK and IWORK arrays, and no error message
143 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
144 *> \endverbatim
145 *>
146 *> \param[out] IWORK
147 *> \verbatim
148 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
149 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
150 *> \endverbatim
151 *>
152 *> \param[in] LIWORK
153 *> \verbatim
154 *> LIWORK is INTEGER
155 *> The dimension of the array IWORK.
156 *> If N <= 1, LIWORK must be at least 1.
157 *> If JOBZ = 'N' and N > 1, LIWORK must be at least 1.
158 *> If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
159 *>
160 *> If LIWORK = -1, then a workspace query is assumed; the
161 *> routine only calculates the optimal sizes of the WORK, RWORK
162 *> and IWORK arrays, returns these values as the first entries
163 *> of the WORK, RWORK and IWORK arrays, and no error message
164 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
165 *> \endverbatim
166 *>
167 *> \param[out] INFO
168 *> \verbatim
169 *> INFO is INTEGER
170 *> = 0: successful exit
171 *> < 0: if INFO = -i, the i-th argument had an illegal value
172 *> > 0: if INFO = i and JOBZ = 'N', then the algorithm failed
173 *> to converge; i off-diagonal elements of an intermediate
174 *> tridiagonal form did not converge to zero;
175 *> if INFO = i and JOBZ = 'V', then the algorithm failed
176 *> to compute an eigenvalue while working on the submatrix
177 *> lying in rows and columns INFO/(N+1) through
178 *> mod(INFO,N+1).
179 *> \endverbatim
180 *
181 * Authors:
182 * ========
183 *
184 *> \author Univ. of Tennessee
185 *> \author Univ. of California Berkeley
186 *> \author Univ. of Colorado Denver
187 *> \author NAG Ltd.
188 *
189 *> \date November 2011
190 *
191 *> \ingroup complex16HEeigen
192 *
193 *> \par Further Details:
194 * =====================
195 *>
196 *> Modified description of INFO. Sven, 16 Feb 05.
197 *
198 *> \par Contributors:
199 * ==================
200 *>
201 *> Jeff Rutter, Computer Science Division, University of California
202 *> at Berkeley, USA
203 *>
204 * =====================================================================
205  SUBROUTINE zheevd( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK,
206  $ lrwork, iwork, liwork, info )
207 *
208 * -- LAPACK driver routine (version 3.4.0) --
209 * -- LAPACK is a software package provided by Univ. of Tennessee, --
210 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
211 * November 2011
212 *
213 * .. Scalar Arguments ..
214  CHARACTER jobz, uplo
215  INTEGER info, lda, liwork, lrwork, lwork, n
216 * ..
217 * .. Array Arguments ..
218  INTEGER iwork( * )
219  DOUBLE PRECISION rwork( * ), w( * )
220  COMPLEX*16 a( lda, * ), work( * )
221 * ..
222 *
223 * =====================================================================
224 *
225 * .. Parameters ..
226  DOUBLE PRECISION zero, one
227  parameter( zero = 0.0d0, one = 1.0d0 )
228  COMPLEX*16 cone
229  parameter( cone = ( 1.0d0, 0.0d0 ) )
230 * ..
231 * .. Local Scalars ..
232  LOGICAL lower, lquery, wantz
233  INTEGER iinfo, imax, inde, indrwk, indtau, indwk2,
234  $ indwrk, iscale, liopt, liwmin, llrwk, llwork,
235  $ llwrk2, lopt, lropt, lrwmin, lwmin
236  DOUBLE PRECISION anrm, bignum, eps, rmax, rmin, safmin, sigma,
237  $ smlnum
238 * ..
239 * .. External Functions ..
240  LOGICAL lsame
241  INTEGER ilaenv
242  DOUBLE PRECISION dlamch, zlanhe
243  EXTERNAL lsame, ilaenv, dlamch, zlanhe
244 * ..
245 * .. External Subroutines ..
246  EXTERNAL dscal, dsterf, xerbla, zhetrd, zlacpy, zlascl,
247  $ zstedc, zunmtr
248 * ..
249 * .. Intrinsic Functions ..
250  INTRINSIC max, sqrt
251 * ..
252 * .. Executable Statements ..
253 *
254 * Test the input parameters.
255 *
256  wantz = lsame( jobz, 'V' )
257  lower = lsame( uplo, 'L' )
258  lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
259 *
260  info = 0
261  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
262  info = -1
263  ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
264  info = -2
265  ELSE IF( n.LT.0 ) THEN
266  info = -3
267  ELSE IF( lda.LT.max( 1, n ) ) THEN
268  info = -5
269  END IF
270 *
271  IF( info.EQ.0 ) THEN
272  IF( n.LE.1 ) THEN
273  lwmin = 1
274  lrwmin = 1
275  liwmin = 1
276  lopt = lwmin
277  lropt = lrwmin
278  liopt = liwmin
279  ELSE
280  IF( wantz ) THEN
281  lwmin = 2*n + n*n
282  lrwmin = 1 + 5*n + 2*n**2
283  liwmin = 3 + 5*n
284  ELSE
285  lwmin = n + 1
286  lrwmin = n
287  liwmin = 1
288  END IF
289  lopt = max( lwmin, n +
290  $ ilaenv( 1, 'ZHETRD', uplo, n, -1, -1, -1 ) )
291  lropt = lrwmin
292  liopt = liwmin
293  END IF
294  work( 1 ) = lopt
295  rwork( 1 ) = lropt
296  iwork( 1 ) = liopt
297 *
298  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
299  info = -8
300  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
301  info = -10
302  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
303  info = -12
304  END IF
305  END IF
306 *
307  IF( info.NE.0 ) THEN
308  CALL xerbla( 'ZHEEVD', -info )
309  return
310  ELSE IF( lquery ) THEN
311  return
312  END IF
313 *
314 * Quick return if possible
315 *
316  IF( n.EQ.0 )
317  $ return
318 *
319  IF( n.EQ.1 ) THEN
320  w( 1 ) = a( 1, 1 )
321  IF( wantz )
322  $ a( 1, 1 ) = cone
323  return
324  END IF
325 *
326 * Get machine constants.
327 *
328  safmin = dlamch( 'Safe minimum' )
329  eps = dlamch( 'Precision' )
330  smlnum = safmin / eps
331  bignum = one / smlnum
332  rmin = sqrt( smlnum )
333  rmax = sqrt( bignum )
334 *
335 * Scale matrix to allowable range, if necessary.
336 *
337  anrm = zlanhe( 'M', uplo, n, a, lda, rwork )
338  iscale = 0
339  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
340  iscale = 1
341  sigma = rmin / anrm
342  ELSE IF( anrm.GT.rmax ) THEN
343  iscale = 1
344  sigma = rmax / anrm
345  END IF
346  IF( iscale.EQ.1 )
347  $ CALL zlascl( uplo, 0, 0, one, sigma, n, n, a, lda, info )
348 *
349 * Call ZHETRD to reduce Hermitian matrix to tridiagonal form.
350 *
351  inde = 1
352  indtau = 1
353  indwrk = indtau + n
354  indrwk = inde + n
355  indwk2 = indwrk + n*n
356  llwork = lwork - indwrk + 1
357  llwrk2 = lwork - indwk2 + 1
358  llrwk = lrwork - indrwk + 1
359  CALL zhetrd( uplo, n, a, lda, w, rwork( inde ), work( indtau ),
360  $ work( indwrk ), llwork, iinfo )
361 *
362 * For eigenvalues only, call DSTERF. For eigenvectors, first call
363 * ZSTEDC to generate the eigenvector matrix, WORK(INDWRK), of the
364 * tridiagonal matrix, then call ZUNMTR to multiply it to the
365 * Householder transformations represented as Householder vectors in
366 * A.
367 *
368  IF( .NOT.wantz ) THEN
369  CALL dsterf( n, w, rwork( inde ), info )
370  ELSE
371  CALL zstedc( 'I', n, w, rwork( inde ), work( indwrk ), n,
372  $ work( indwk2 ), llwrk2, rwork( indrwk ), llrwk,
373  $ iwork, liwork, info )
374  CALL zunmtr( 'L', uplo, 'N', n, n, a, lda, work( indtau ),
375  $ work( indwrk ), n, work( indwk2 ), llwrk2, iinfo )
376  CALL zlacpy( 'A', n, n, work( indwrk ), n, a, lda )
377  END IF
378 *
379 * If matrix was scaled, then rescale eigenvalues appropriately.
380 *
381  IF( iscale.EQ.1 ) THEN
382  IF( info.EQ.0 ) THEN
383  imax = n
384  ELSE
385  imax = info - 1
386  END IF
387  CALL dscal( imax, one / sigma, w, 1 )
388  END IF
389 *
390  work( 1 ) = lopt
391  rwork( 1 ) = lropt
392  iwork( 1 ) = liopt
393 *
394  return
395 *
396 * End of ZHEEVD
397 *
398  END