LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zhpevx.f
Go to the documentation of this file.
1 *> \brief <b> ZHPEVX 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 ZHPEVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhpevx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhpevx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhpevx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZHPEVX( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU,
22 * ABSTOL, M, W, Z, LDZ, WORK, RWORK, IWORK,
23 * IFAIL, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBZ, RANGE, UPLO
27 * INTEGER IL, INFO, IU, LDZ, M, N
28 * DOUBLE PRECISION ABSTOL, VL, VU
29 * ..
30 * .. Array Arguments ..
31 * INTEGER IFAIL( * ), IWORK( * )
32 * DOUBLE PRECISION RWORK( * ), W( * )
33 * COMPLEX*16 AP( * ), WORK( * ), Z( LDZ, * )
34 * ..
35 *
36 *
37 *> \par Purpose:
38 * =============
39 *>
40 *> \verbatim
41 *>
42 *> ZHPEVX computes selected eigenvalues and, optionally, eigenvectors
43 *> of a complex Hermitian matrix A in packed storage.
44 *> Eigenvalues/vectors can be selected by specifying either a range of
45 *> values or a range of indices for the desired eigenvalues.
46 *> \endverbatim
47 *
48 * Arguments:
49 * ==========
50 *
51 *> \param[in] JOBZ
52 *> \verbatim
53 *> JOBZ is CHARACTER*1
54 *> = 'N': Compute eigenvalues only;
55 *> = 'V': Compute eigenvalues and eigenvectors.
56 *> \endverbatim
57 *>
58 *> \param[in] RANGE
59 *> \verbatim
60 *> RANGE is CHARACTER*1
61 *> = 'A': all eigenvalues will be found;
62 *> = 'V': all eigenvalues in the half-open interval (VL,VU]
63 *> will be found;
64 *> = 'I': the IL-th through IU-th eigenvalues will be found.
65 *> \endverbatim
66 *>
67 *> \param[in] UPLO
68 *> \verbatim
69 *> UPLO is CHARACTER*1
70 *> = 'U': Upper triangle of A is stored;
71 *> = 'L': Lower triangle of A is stored.
72 *> \endverbatim
73 *>
74 *> \param[in] N
75 *> \verbatim
76 *> N is INTEGER
77 *> The order of the matrix A. N >= 0.
78 *> \endverbatim
79 *>
80 *> \param[in,out] AP
81 *> \verbatim
82 *> AP is COMPLEX*16 array, dimension (N*(N+1)/2)
83 *> On entry, the upper or lower triangle of the Hermitian matrix
84 *> A, packed columnwise in a linear array. The j-th column of A
85 *> is stored in the array AP as follows:
86 *> if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
87 *> if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
88 *>
89 *> On exit, AP is overwritten by values generated during the
90 *> reduction to tridiagonal form. If UPLO = 'U', the diagonal
91 *> and first superdiagonal of the tridiagonal matrix T overwrite
92 *> the corresponding elements of A, and if UPLO = 'L', the
93 *> diagonal and first subdiagonal of T overwrite the
94 *> corresponding elements of A.
95 *> \endverbatim
96 *>
97 *> \param[in] VL
98 *> \verbatim
99 *> VL is DOUBLE PRECISION
100 *> \endverbatim
101 *>
102 *> \param[in] VU
103 *> \verbatim
104 *> VU is DOUBLE PRECISION
105 *> If RANGE='V', the lower and upper bounds of the interval to
106 *> be searched for eigenvalues. VL < VU.
107 *> Not referenced if RANGE = 'A' or 'I'.
108 *> \endverbatim
109 *>
110 *> \param[in] IL
111 *> \verbatim
112 *> IL is INTEGER
113 *> \endverbatim
114 *>
115 *> \param[in] IU
116 *> \verbatim
117 *> IU is INTEGER
118 *> If RANGE='I', the indices (in ascending order) of the
119 *> smallest and largest eigenvalues to be returned.
120 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
121 *> Not referenced if RANGE = 'A' or 'V'.
122 *> \endverbatim
123 *>
124 *> \param[in] ABSTOL
125 *> \verbatim
126 *> ABSTOL is DOUBLE PRECISION
127 *> The absolute error tolerance for the eigenvalues.
128 *> An approximate eigenvalue is accepted as converged
129 *> when it is determined to lie in an interval [a,b]
130 *> of width less than or equal to
131 *>
132 *> ABSTOL + EPS * max( |a|,|b| ) ,
133 *>
134 *> where EPS is the machine precision. If ABSTOL is less than
135 *> or equal to zero, then EPS*|T| will be used in its place,
136 *> where |T| is the 1-norm of the tridiagonal matrix obtained
137 *> by reducing AP to tridiagonal form.
138 *>
139 *> Eigenvalues will be computed most accurately when ABSTOL is
140 *> set to twice the underflow threshold 2*DLAMCH('S'), not zero.
141 *> If this routine returns with INFO>0, indicating that some
142 *> eigenvectors did not converge, try setting ABSTOL to
143 *> 2*DLAMCH('S').
144 *>
145 *> See "Computing Small Singular Values of Bidiagonal Matrices
146 *> with Guaranteed High Relative Accuracy," by Demmel and
147 *> Kahan, LAPACK Working Note #3.
148 *> \endverbatim
149 *>
150 *> \param[out] M
151 *> \verbatim
152 *> M is INTEGER
153 *> The total number of eigenvalues found. 0 <= M <= N.
154 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
155 *> \endverbatim
156 *>
157 *> \param[out] W
158 *> \verbatim
159 *> W is DOUBLE PRECISION array, dimension (N)
160 *> If INFO = 0, the selected eigenvalues in ascending order.
161 *> \endverbatim
162 *>
163 *> \param[out] Z
164 *> \verbatim
165 *> Z is COMPLEX*16 array, dimension (LDZ, max(1,M))
166 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
167 *> contain the orthonormal eigenvectors of the matrix A
168 *> corresponding to the selected eigenvalues, with the i-th
169 *> column of Z holding the eigenvector associated with W(i).
170 *> If an eigenvector fails to converge, then that column of Z
171 *> contains the latest approximation to the eigenvector, and
172 *> the index of the eigenvector is returned in IFAIL.
173 *> If JOBZ = 'N', then Z is not referenced.
174 *> Note: the user must ensure that at least max(1,M) columns are
175 *> supplied in the array Z; if RANGE = 'V', the exact value of M
176 *> is not known in advance and an upper bound must be used.
177 *> \endverbatim
178 *>
179 *> \param[in] LDZ
180 *> \verbatim
181 *> LDZ is INTEGER
182 *> The leading dimension of the array Z. LDZ >= 1, and if
183 *> JOBZ = 'V', LDZ >= max(1,N).
184 *> \endverbatim
185 *>
186 *> \param[out] WORK
187 *> \verbatim
188 *> WORK is COMPLEX*16 array, dimension (2*N)
189 *> \endverbatim
190 *>
191 *> \param[out] RWORK
192 *> \verbatim
193 *> RWORK is DOUBLE PRECISION array, dimension (7*N)
194 *> \endverbatim
195 *>
196 *> \param[out] IWORK
197 *> \verbatim
198 *> IWORK is INTEGER array, dimension (5*N)
199 *> \endverbatim
200 *>
201 *> \param[out] IFAIL
202 *> \verbatim
203 *> IFAIL is INTEGER array, dimension (N)
204 *> If JOBZ = 'V', then if INFO = 0, the first M elements of
205 *> IFAIL are zero. If INFO > 0, then IFAIL contains the
206 *> indices of the eigenvectors that failed to converge.
207 *> If JOBZ = 'N', then IFAIL is not referenced.
208 *> \endverbatim
209 *>
210 *> \param[out] INFO
211 *> \verbatim
212 *> INFO is INTEGER
213 *> = 0: successful exit
214 *> < 0: if INFO = -i, the i-th argument had an illegal value
215 *> > 0: if INFO = i, then i eigenvectors failed to converge.
216 *> Their indices are stored in array IFAIL.
217 *> \endverbatim
218 *
219 * Authors:
220 * ========
221 *
222 *> \author Univ. of Tennessee
223 *> \author Univ. of California Berkeley
224 *> \author Univ. of Colorado Denver
225 *> \author NAG Ltd.
226 *
227 *> \date November 2011
228 *
229 *> \ingroup complex16OTHEReigen
230 *
231 * =====================================================================
232  SUBROUTINE zhpevx( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU,
233  $ abstol, m, w, z, ldz, work, rwork, iwork,
234  $ ifail, info )
235 *
236 * -- LAPACK driver routine (version 3.4.0) --
237 * -- LAPACK is a software package provided by Univ. of Tennessee, --
238 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
239 * November 2011
240 *
241 * .. Scalar Arguments ..
242  CHARACTER jobz, range, uplo
243  INTEGER il, info, iu, ldz, m, n
244  DOUBLE PRECISION abstol, vl, vu
245 * ..
246 * .. Array Arguments ..
247  INTEGER ifail( * ), iwork( * )
248  DOUBLE PRECISION rwork( * ), w( * )
249  COMPLEX*16 ap( * ), work( * ), z( ldz, * )
250 * ..
251 *
252 * =====================================================================
253 *
254 * .. Parameters ..
255  DOUBLE PRECISION zero, one
256  parameter( zero = 0.0d0, one = 1.0d0 )
257  COMPLEX*16 cone
258  parameter( cone = ( 1.0d0, 0.0d0 ) )
259 * ..
260 * .. Local Scalars ..
261  LOGICAL alleig, indeig, test, valeig, wantz
262  CHARACTER order
263  INTEGER i, iinfo, imax, indd, inde, indee, indibl,
264  $ indisp, indiwk, indrwk, indtau, indwrk, iscale,
265  $ itmp1, j, jj, nsplit
266  DOUBLE PRECISION abstll, anrm, bignum, eps, rmax, rmin, safmin,
267  $ sigma, smlnum, tmp1, vll, vuu
268 * ..
269 * .. External Functions ..
270  LOGICAL lsame
271  DOUBLE PRECISION dlamch, zlanhp
272  EXTERNAL lsame, dlamch, zlanhp
273 * ..
274 * .. External Subroutines ..
275  EXTERNAL dcopy, dscal, dstebz, dsterf, xerbla, zdscal,
276  $ zhptrd, zstein, zsteqr, zswap, zupgtr, zupmtr
277 * ..
278 * .. Intrinsic Functions ..
279  INTRINSIC dble, max, min, sqrt
280 * ..
281 * .. Executable Statements ..
282 *
283 * Test the input parameters.
284 *
285  wantz = lsame( jobz, 'V' )
286  alleig = lsame( range, 'A' )
287  valeig = lsame( range, 'V' )
288  indeig = lsame( range, 'I' )
289 *
290  info = 0
291  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
292  info = -1
293  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
294  info = -2
295  ELSE IF( .NOT.( lsame( uplo, 'L' ) .OR. lsame( uplo, 'U' ) ) )
296  $ THEN
297  info = -3
298  ELSE IF( n.LT.0 ) THEN
299  info = -4
300  ELSE
301  IF( valeig ) THEN
302  IF( n.GT.0 .AND. vu.LE.vl )
303  $ info = -7
304  ELSE IF( indeig ) THEN
305  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
306  info = -8
307  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
308  info = -9
309  END IF
310  END IF
311  END IF
312  IF( info.EQ.0 ) THEN
313  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
314  $ info = -14
315  END IF
316 *
317  IF( info.NE.0 ) THEN
318  CALL xerbla( 'ZHPEVX', -info )
319  return
320  END IF
321 *
322 * Quick return if possible
323 *
324  m = 0
325  IF( n.EQ.0 )
326  $ return
327 *
328  IF( n.EQ.1 ) THEN
329  IF( alleig .OR. indeig ) THEN
330  m = 1
331  w( 1 ) = ap( 1 )
332  ELSE
333  IF( vl.LT.dble( ap( 1 ) ) .AND. vu.GE.dble( ap( 1 ) ) ) THEN
334  m = 1
335  w( 1 ) = ap( 1 )
336  END IF
337  END IF
338  IF( wantz )
339  $ z( 1, 1 ) = cone
340  return
341  END IF
342 *
343 * Get machine constants.
344 *
345  safmin = dlamch( 'Safe minimum' )
346  eps = dlamch( 'Precision' )
347  smlnum = safmin / eps
348  bignum = one / smlnum
349  rmin = sqrt( smlnum )
350  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
351 *
352 * Scale matrix to allowable range, if necessary.
353 *
354  iscale = 0
355  abstll = abstol
356  IF( valeig ) THEN
357  vll = vl
358  vuu = vu
359  ELSE
360  vll = zero
361  vuu = zero
362  END IF
363  anrm = zlanhp( 'M', uplo, n, ap, rwork )
364  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
365  iscale = 1
366  sigma = rmin / anrm
367  ELSE IF( anrm.GT.rmax ) THEN
368  iscale = 1
369  sigma = rmax / anrm
370  END IF
371  IF( iscale.EQ.1 ) THEN
372  CALL zdscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
373  IF( abstol.GT.0 )
374  $ abstll = abstol*sigma
375  IF( valeig ) THEN
376  vll = vl*sigma
377  vuu = vu*sigma
378  END IF
379  END IF
380 *
381 * Call ZHPTRD to reduce Hermitian packed matrix to tridiagonal form.
382 *
383  indd = 1
384  inde = indd + n
385  indrwk = inde + n
386  indtau = 1
387  indwrk = indtau + n
388  CALL zhptrd( uplo, n, ap, rwork( indd ), rwork( inde ),
389  $ work( indtau ), iinfo )
390 *
391 * If all eigenvalues are desired and ABSTOL is less than or equal
392 * to zero, then call DSTERF or ZUPGTR and ZSTEQR. If this fails
393 * for some eigenvalue, then try DSTEBZ.
394 *
395  test = .false.
396  IF (indeig) THEN
397  IF (il.EQ.1 .AND. iu.EQ.n) THEN
398  test = .true.
399  END IF
400  END IF
401  IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
402  CALL dcopy( n, rwork( indd ), 1, w, 1 )
403  indee = indrwk + 2*n
404  IF( .NOT.wantz ) THEN
405  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
406  CALL dsterf( n, w, rwork( indee ), info )
407  ELSE
408  CALL zupgtr( uplo, n, ap, work( indtau ), z, ldz,
409  $ work( indwrk ), iinfo )
410  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
411  CALL zsteqr( jobz, n, w, rwork( indee ), z, ldz,
412  $ rwork( indrwk ), info )
413  IF( info.EQ.0 ) THEN
414  DO 10 i = 1, n
415  ifail( i ) = 0
416  10 continue
417  END IF
418  END IF
419  IF( info.EQ.0 ) THEN
420  m = n
421  go to 20
422  END IF
423  info = 0
424  END IF
425 *
426 * Otherwise, call DSTEBZ and, if eigenvectors are desired, ZSTEIN.
427 *
428  IF( wantz ) THEN
429  order = 'B'
430  ELSE
431  order = 'E'
432  END IF
433  indibl = 1
434  indisp = indibl + n
435  indiwk = indisp + n
436  CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
437  $ rwork( indd ), rwork( inde ), m, nsplit, w,
438  $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
439  $ iwork( indiwk ), info )
440 *
441  IF( wantz ) THEN
442  CALL zstein( n, rwork( indd ), rwork( inde ), m, w,
443  $ iwork( indibl ), iwork( indisp ), z, ldz,
444  $ rwork( indrwk ), iwork( indiwk ), ifail, info )
445 *
446 * Apply unitary matrix used in reduction to tridiagonal
447 * form to eigenvectors returned by ZSTEIN.
448 *
449  indwrk = indtau + n
450  CALL zupmtr( 'L', uplo, 'N', n, m, ap, work( indtau ), z, ldz,
451  $ work( indwrk ), iinfo )
452  END IF
453 *
454 * If matrix was scaled, then rescale eigenvalues appropriately.
455 *
456  20 continue
457  IF( iscale.EQ.1 ) THEN
458  IF( info.EQ.0 ) THEN
459  imax = m
460  ELSE
461  imax = info - 1
462  END IF
463  CALL dscal( imax, one / sigma, w, 1 )
464  END IF
465 *
466 * If eigenvalues are not in order, then sort them, along with
467 * eigenvectors.
468 *
469  IF( wantz ) THEN
470  DO 40 j = 1, m - 1
471  i = 0
472  tmp1 = w( j )
473  DO 30 jj = j + 1, m
474  IF( w( jj ).LT.tmp1 ) THEN
475  i = jj
476  tmp1 = w( jj )
477  END IF
478  30 continue
479 *
480  IF( i.NE.0 ) THEN
481  itmp1 = iwork( indibl+i-1 )
482  w( i ) = w( j )
483  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
484  w( j ) = tmp1
485  iwork( indibl+j-1 ) = itmp1
486  CALL zswap( n, z( 1, i ), 1, z( 1, j ), 1 )
487  IF( info.NE.0 ) THEN
488  itmp1 = ifail( i )
489  ifail( i ) = ifail( j )
490  ifail( j ) = itmp1
491  END IF
492  END IF
493  40 continue
494  END IF
495 *
496  return
497 *
498 * End of ZHPEVX
499 *
500  END