LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zhpevd.f
Go to the documentation of this file.
1 *> \brief <b> ZHPEVD computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER 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 ZHPEVD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhpevd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhpevd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhpevd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHPEVD( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK,
22 * RWORK, LRWORK, IWORK, LIWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER JOBZ, UPLO
26 * INTEGER INFO, LDZ, LIWORK, LRWORK, LWORK, N
27 * ..
28 * .. Array Arguments ..
29 * INTEGER IWORK( * )
30 * DOUBLE PRECISION RWORK( * ), W( * )
31 * COMPLEX*16 AP( * ), WORK( * ), Z( LDZ, * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> ZHPEVD computes all the eigenvalues and, optionally, eigenvectors of
41 *> a complex Hermitian matrix A in packed storage. If eigenvectors are
42 *> desired, it uses a 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] AP
76 *> \verbatim
77 *> AP is COMPLEX*16 array, dimension (N*(N+1)/2)
78 *> On entry, the upper or lower triangle of the Hermitian matrix
79 *> A, packed columnwise in a linear array. The j-th column of A
80 *> is stored in the array AP as follows:
81 *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
82 *> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
83 *>
84 *> On exit, AP is overwritten by values generated during the
85 *> reduction to tridiagonal form. If UPLO = 'U', the diagonal
86 *> and first superdiagonal of the tridiagonal matrix T overwrite
87 *> the corresponding elements of A, and if UPLO = 'L', the
88 *> diagonal and first subdiagonal of T overwrite the
89 *> corresponding elements of A.
90 *> \endverbatim
91 *>
92 *> \param[out] W
93 *> \verbatim
94 *> W is DOUBLE PRECISION array, dimension (N)
95 *> If INFO = 0, the eigenvalues in ascending order.
96 *> \endverbatim
97 *>
98 *> \param[out] Z
99 *> \verbatim
100 *> Z is COMPLEX*16 array, dimension (LDZ, N)
101 *> If JOBZ = 'V', then if INFO = 0, Z contains the orthonormal
102 *> eigenvectors of the matrix A, with the i-th column of Z
103 *> holding the eigenvector associated with W(i).
104 *> If JOBZ = 'N', then Z is not referenced.
105 *> \endverbatim
106 *>
107 *> \param[in] LDZ
108 *> \verbatim
109 *> LDZ is INTEGER
110 *> The leading dimension of the array Z. LDZ >= 1, and if
111 *> JOBZ = 'V', LDZ >= max(1,N).
112 *> \endverbatim
113 *>
114 *> \param[out] WORK
115 *> \verbatim
116 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
117 *> On exit, if INFO = 0, WORK(1) returns the required LWORK.
118 *> \endverbatim
119 *>
120 *> \param[in] LWORK
121 *> \verbatim
122 *> LWORK is INTEGER
123 *> The dimension of array WORK.
124 *> If N <= 1, LWORK must be at least 1.
125 *> If JOBZ = 'N' and N > 1, LWORK must be at least N.
126 *> If JOBZ = 'V' and N > 1, LWORK must be at least 2*N.
127 *>
128 *> If LWORK = -1, then a workspace query is assumed; the routine
129 *> only calculates the required sizes of the WORK, RWORK and
130 *> IWORK arrays, returns these values as the first entries of
131 *> the WORK, RWORK and IWORK arrays, and no error message
132 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
133 *> \endverbatim
134 *>
135 *> \param[out] RWORK
136 *> \verbatim
137 *> RWORK is DOUBLE PRECISION array,
138 *> dimension (LRWORK)
139 *> On exit, if INFO = 0, RWORK(1) returns the required LRWORK.
140 *> \endverbatim
141 *>
142 *> \param[in] LRWORK
143 *> \verbatim
144 *> LRWORK is INTEGER
145 *> The dimension of array RWORK.
146 *> If N <= 1, LRWORK must be at least 1.
147 *> If JOBZ = 'N' and N > 1, LRWORK must be at least N.
148 *> If JOBZ = 'V' and N > 1, LRWORK must be at least
149 *> 1 + 5*N + 2*N**2.
150 *>
151 *> If LRWORK = -1, then a workspace query is assumed; the
152 *> routine only calculates the required sizes of the WORK, RWORK
153 *> and IWORK arrays, returns these values as the first entries
154 *> of the WORK, RWORK and IWORK arrays, and no error message
155 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
156 *> \endverbatim
157 *>
158 *> \param[out] IWORK
159 *> \verbatim
160 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
161 *> On exit, if INFO = 0, IWORK(1) returns the required LIWORK.
162 *> \endverbatim
163 *>
164 *> \param[in] LIWORK
165 *> \verbatim
166 *> LIWORK is INTEGER
167 *> The dimension of array IWORK.
168 *> If JOBZ = 'N' or N <= 1, LIWORK must be at least 1.
169 *> If JOBZ = 'V' and N > 1, LIWORK must be at least 3 + 5*N.
170 *>
171 *> If LIWORK = -1, then a workspace query is assumed; the
172 *> routine only calculates the required sizes of the WORK, RWORK
173 *> and IWORK arrays, returns these values as the first entries
174 *> of the WORK, RWORK and IWORK arrays, and no error message
175 *> related to LWORK or LRWORK or LIWORK is issued by XERBLA.
176 *> \endverbatim
177 *>
178 *> \param[out] INFO
179 *> \verbatim
180 *> INFO is INTEGER
181 *> = 0: successful exit
182 *> < 0: if INFO = -i, the i-th argument had an illegal value.
183 *> > 0: if INFO = i, the algorithm failed to converge; i
184 *> off-diagonal elements of an intermediate tridiagonal
185 *> form did not converge to zero.
186 *> \endverbatim
187 *
188 * Authors:
189 * ========
190 *
191 *> \author Univ. of Tennessee
192 *> \author Univ. of California Berkeley
193 *> \author Univ. of Colorado Denver
194 *> \author NAG Ltd.
195 *
196 *> \date November 2011
197 *
198 *> \ingroup complex16OTHEReigen
199 *
200 * =====================================================================
201  SUBROUTINE zhpevd( JOBZ, UPLO, N, AP, W, Z, LDZ, WORK, LWORK,
202  $ rwork, lrwork, iwork, liwork, info )
203 *
204 * -- LAPACK driver routine (version 3.4.0) --
205 * -- LAPACK is a software package provided by Univ. of Tennessee, --
206 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
207 * November 2011
208 *
209 * .. Scalar Arguments ..
210  CHARACTER jobz, uplo
211  INTEGER info, ldz, liwork, lrwork, lwork, n
212 * ..
213 * .. Array Arguments ..
214  INTEGER iwork( * )
215  DOUBLE PRECISION rwork( * ), w( * )
216  COMPLEX*16 ap( * ), work( * ), z( ldz, * )
217 * ..
218 *
219 * =====================================================================
220 *
221 * .. Parameters ..
222  DOUBLE PRECISION zero, one
223  parameter( zero = 0.0d+0, one = 1.0d+0 )
224  COMPLEX*16 cone
225  parameter( cone = ( 1.0d+0, 0.0d+0 ) )
226 * ..
227 * .. Local Scalars ..
228  LOGICAL lquery, wantz
229  INTEGER iinfo, imax, inde, indrwk, indtau, indwrk,
230  $ iscale, liwmin, llrwk, llwrk, lrwmin, lwmin
231  DOUBLE PRECISION anrm, bignum, eps, rmax, rmin, safmin, sigma,
232  $ smlnum
233 * ..
234 * .. External Functions ..
235  LOGICAL lsame
236  DOUBLE PRECISION dlamch, zlanhp
237  EXTERNAL lsame, dlamch, zlanhp
238 * ..
239 * .. External Subroutines ..
240  EXTERNAL dscal, dsterf, xerbla, zdscal, zhptrd, zstedc,
241  $ zupmtr
242 * ..
243 * .. Intrinsic Functions ..
244  INTRINSIC sqrt
245 * ..
246 * .. Executable Statements ..
247 *
248 * Test the input parameters.
249 *
250  wantz = lsame( jobz, 'V' )
251  lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 .OR. liwork.EQ.-1 )
252 *
253  info = 0
254  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
255  info = -1
256  ELSE IF( .NOT.( lsame( uplo, 'L' ) .OR. lsame( uplo, 'U' ) ) )
257  $ THEN
258  info = -2
259  ELSE IF( n.LT.0 ) THEN
260  info = -3
261  ELSE IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
262  info = -7
263  END IF
264 *
265  IF( info.EQ.0 ) THEN
266  IF( n.LE.1 ) THEN
267  lwmin = 1
268  liwmin = 1
269  lrwmin = 1
270  ELSE
271  IF( wantz ) THEN
272  lwmin = 2*n
273  lrwmin = 1 + 5*n + 2*n**2
274  liwmin = 3 + 5*n
275  ELSE
276  lwmin = n
277  lrwmin = n
278  liwmin = 1
279  END IF
280  END IF
281  work( 1 ) = lwmin
282  rwork( 1 ) = lrwmin
283  iwork( 1 ) = liwmin
284 *
285  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
286  info = -9
287  ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
288  info = -11
289  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
290  info = -13
291  END IF
292  END IF
293 *
294  IF( info.NE.0 ) THEN
295  CALL xerbla( 'ZHPEVD', -info )
296  return
297  ELSE IF( lquery ) THEN
298  return
299  END IF
300 *
301 * Quick return if possible
302 *
303  IF( n.EQ.0 )
304  $ return
305 *
306  IF( n.EQ.1 ) THEN
307  w( 1 ) = ap( 1 )
308  IF( wantz )
309  $ z( 1, 1 ) = cone
310  return
311  END IF
312 *
313 * Get machine constants.
314 *
315  safmin = dlamch( 'Safe minimum' )
316  eps = dlamch( 'Precision' )
317  smlnum = safmin / eps
318  bignum = one / smlnum
319  rmin = sqrt( smlnum )
320  rmax = sqrt( bignum )
321 *
322 * Scale matrix to allowable range, if necessary.
323 *
324  anrm = zlanhp( 'M', uplo, n, ap, rwork )
325  iscale = 0
326  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
327  iscale = 1
328  sigma = rmin / anrm
329  ELSE IF( anrm.GT.rmax ) THEN
330  iscale = 1
331  sigma = rmax / anrm
332  END IF
333  IF( iscale.EQ.1 ) THEN
334  CALL zdscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
335  END IF
336 *
337 * Call ZHPTRD to reduce Hermitian packed matrix to tridiagonal form.
338 *
339  inde = 1
340  indtau = 1
341  indrwk = inde + n
342  indwrk = indtau + n
343  llwrk = lwork - indwrk + 1
344  llrwk = lrwork - indrwk + 1
345  CALL zhptrd( uplo, n, ap, w, rwork( inde ), work( indtau ),
346  $ iinfo )
347 *
348 * For eigenvalues only, call DSTERF. For eigenvectors, first call
349 * ZUPGTR to generate the orthogonal matrix, then call ZSTEDC.
350 *
351  IF( .NOT.wantz ) THEN
352  CALL dsterf( n, w, rwork( inde ), info )
353  ELSE
354  CALL zstedc( 'I', n, w, rwork( inde ), z, ldz, work( indwrk ),
355  $ llwrk, rwork( indrwk ), llrwk, iwork, liwork,
356  $ info )
357  CALL zupmtr( 'L', uplo, 'N', n, n, ap, work( indtau ), z, ldz,
358  $ work( indwrk ), iinfo )
359  END IF
360 *
361 * If matrix was scaled, then rescale eigenvalues appropriately.
362 *
363  IF( iscale.EQ.1 ) THEN
364  IF( info.EQ.0 ) THEN
365  imax = n
366  ELSE
367  imax = info - 1
368  END IF
369  CALL dscal( imax, one / sigma, w, 1 )
370  END IF
371 *
372  work( 1 ) = lwmin
373  rwork( 1 ) = lrwmin
374  iwork( 1 ) = liwmin
375  return
376 *
377 * End of ZHPEVD
378 *
379  END