LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 *> If RANGE='V', the lower bound of the interval to
101 *> be searched for eigenvalues. VL < VU.
102 *> Not referenced if RANGE = 'A' or 'I'.
103 *> \endverbatim
104 *>
105 *> \param[in] VU
106 *> \verbatim
107 *> VU is DOUBLE PRECISION
108 *> If RANGE='V', the upper bound of the interval to
109 *> be searched for eigenvalues. VL < VU.
110 *> Not referenced if RANGE = 'A' or 'I'.
111 *> \endverbatim
112 *>
113 *> \param[in] IL
114 *> \verbatim
115 *> IL is INTEGER
116 *> If RANGE='I', the index of the
117 *> smallest eigenvalue to be returned.
118 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
119 *> Not referenced if RANGE = 'A' or 'V'.
120 *> \endverbatim
121 *>
122 *> \param[in] IU
123 *> \verbatim
124 *> IU is INTEGER
125 *> If RANGE='I', the index of the
126 *> largest eigenvalue to be returned.
127 *> 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
128 *> Not referenced if RANGE = 'A' or 'V'.
129 *> \endverbatim
130 *>
131 *> \param[in] ABSTOL
132 *> \verbatim
133 *> ABSTOL is DOUBLE PRECISION
134 *> The absolute error tolerance for the eigenvalues.
135 *> An approximate eigenvalue is accepted as converged
136 *> when it is determined to lie in an interval [a,b]
137 *> of width less than or equal to
138 *>
139 *> ABSTOL + EPS * max( |a|,|b| ) ,
140 *>
141 *> where EPS is the machine precision. If ABSTOL is less than
142 *> or equal to zero, then EPS*|T| will be used in its place,
143 *> where |T| is the 1-norm of the tridiagonal matrix obtained
144 *> by reducing AP to tridiagonal form.
145 *>
146 *> Eigenvalues will be computed most accurately when ABSTOL is
147 *> set to twice the underflow threshold 2*DLAMCH('S'), not zero.
148 *> If this routine returns with INFO>0, indicating that some
149 *> eigenvectors did not converge, try setting ABSTOL to
150 *> 2*DLAMCH('S').
151 *>
152 *> See "Computing Small Singular Values of Bidiagonal Matrices
153 *> with Guaranteed High Relative Accuracy," by Demmel and
154 *> Kahan, LAPACK Working Note #3.
155 *> \endverbatim
156 *>
157 *> \param[out] M
158 *> \verbatim
159 *> M is INTEGER
160 *> The total number of eigenvalues found. 0 <= M <= N.
161 *> If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
162 *> \endverbatim
163 *>
164 *> \param[out] W
165 *> \verbatim
166 *> W is DOUBLE PRECISION array, dimension (N)
167 *> If INFO = 0, the selected eigenvalues in ascending order.
168 *> \endverbatim
169 *>
170 *> \param[out] Z
171 *> \verbatim
172 *> Z is COMPLEX*16 array, dimension (LDZ, max(1,M))
173 *> If JOBZ = 'V', then if INFO = 0, the first M columns of Z
174 *> contain the orthonormal eigenvectors of the matrix A
175 *> corresponding to the selected eigenvalues, with the i-th
176 *> column of Z holding the eigenvector associated with W(i).
177 *> If an eigenvector fails to converge, then that column of Z
178 *> contains the latest approximation to the eigenvector, and
179 *> the index of the eigenvector is returned in IFAIL.
180 *> If JOBZ = 'N', then Z is not referenced.
181 *> Note: the user must ensure that at least max(1,M) columns are
182 *> supplied in the array Z; if RANGE = 'V', the exact value of M
183 *> is not known in advance and an upper bound must be used.
184 *> \endverbatim
185 *>
186 *> \param[in] LDZ
187 *> \verbatim
188 *> LDZ is INTEGER
189 *> The leading dimension of the array Z. LDZ >= 1, and if
190 *> JOBZ = 'V', LDZ >= max(1,N).
191 *> \endverbatim
192 *>
193 *> \param[out] WORK
194 *> \verbatim
195 *> WORK is COMPLEX*16 array, dimension (2*N)
196 *> \endverbatim
197 *>
198 *> \param[out] RWORK
199 *> \verbatim
200 *> RWORK is DOUBLE PRECISION array, dimension (7*N)
201 *> \endverbatim
202 *>
203 *> \param[out] IWORK
204 *> \verbatim
205 *> IWORK is INTEGER array, dimension (5*N)
206 *> \endverbatim
207 *>
208 *> \param[out] IFAIL
209 *> \verbatim
210 *> IFAIL is INTEGER array, dimension (N)
211 *> If JOBZ = 'V', then if INFO = 0, the first M elements of
212 *> IFAIL are zero. If INFO > 0, then IFAIL contains the
213 *> indices of the eigenvectors that failed to converge.
214 *> If JOBZ = 'N', then IFAIL is not referenced.
215 *> \endverbatim
216 *>
217 *> \param[out] INFO
218 *> \verbatim
219 *> INFO is INTEGER
220 *> = 0: successful exit
221 *> < 0: if INFO = -i, the i-th argument had an illegal value
222 *> > 0: if INFO = i, then i eigenvectors failed to converge.
223 *> Their indices are stored in array IFAIL.
224 *> \endverbatim
225 *
226 * Authors:
227 * ========
228 *
229 *> \author Univ. of Tennessee
230 *> \author Univ. of California Berkeley
231 *> \author Univ. of Colorado Denver
232 *> \author NAG Ltd.
233 *
234 *> \ingroup complex16OTHEReigen
235 *
236 * =====================================================================
237  SUBROUTINE zhpevx( JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU,
238  $ ABSTOL, M, W, Z, LDZ, WORK, RWORK, IWORK,
239  $ IFAIL, INFO )
240 *
241 * -- LAPACK driver routine --
242 * -- LAPACK is a software package provided by Univ. of Tennessee, --
243 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
244 *
245 * .. Scalar Arguments ..
246  CHARACTER JOBZ, RANGE, UPLO
247  INTEGER IL, INFO, IU, LDZ, M, N
248  DOUBLE PRECISION ABSTOL, VL, VU
249 * ..
250 * .. Array Arguments ..
251  INTEGER IFAIL( * ), IWORK( * )
252  DOUBLE PRECISION RWORK( * ), W( * )
253  COMPLEX*16 AP( * ), WORK( * ), Z( LDZ, * )
254 * ..
255 *
256 * =====================================================================
257 *
258 * .. Parameters ..
259  DOUBLE PRECISION ZERO, ONE
260  PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
261  COMPLEX*16 CONE
262  parameter( cone = ( 1.0d0, 0.0d0 ) )
263 * ..
264 * .. Local Scalars ..
265  LOGICAL ALLEIG, INDEIG, TEST, VALEIG, WANTZ
266  CHARACTER ORDER
267  INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
268  $ indisp, indiwk, indrwk, indtau, indwrk, iscale,
269  $ itmp1, j, jj, nsplit
270  DOUBLE PRECISION ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
271  $ SIGMA, SMLNUM, TMP1, VLL, VUU
272 * ..
273 * .. External Functions ..
274  LOGICAL LSAME
275  DOUBLE PRECISION DLAMCH, ZLANHP
276  EXTERNAL lsame, dlamch, zlanhp
277 * ..
278 * .. External Subroutines ..
279  EXTERNAL dcopy, dscal, dstebz, dsterf, xerbla, zdscal,
281 * ..
282 * .. Intrinsic Functions ..
283  INTRINSIC dble, max, min, sqrt
284 * ..
285 * .. Executable Statements ..
286 *
287 * Test the input parameters.
288 *
289  wantz = lsame( jobz, 'V' )
290  alleig = lsame( range, 'A' )
291  valeig = lsame( range, 'V' )
292  indeig = lsame( range, 'I' )
293 *
294  info = 0
295  IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
296  info = -1
297  ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
298  info = -2
299  ELSE IF( .NOT.( lsame( uplo, 'L' ) .OR. lsame( uplo, 'U' ) ) )
300  $ THEN
301  info = -3
302  ELSE IF( n.LT.0 ) THEN
303  info = -4
304  ELSE
305  IF( valeig ) THEN
306  IF( n.GT.0 .AND. vu.LE.vl )
307  $ info = -7
308  ELSE IF( indeig ) THEN
309  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
310  info = -8
311  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
312  info = -9
313  END IF
314  END IF
315  END IF
316  IF( info.EQ.0 ) THEN
317  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
318  $ info = -14
319  END IF
320 *
321  IF( info.NE.0 ) THEN
322  CALL xerbla( 'ZHPEVX', -info )
323  RETURN
324  END IF
325 *
326 * Quick return if possible
327 *
328  m = 0
329  IF( n.EQ.0 )
330  $ RETURN
331 *
332  IF( n.EQ.1 ) THEN
333  IF( alleig .OR. indeig ) THEN
334  m = 1
335  w( 1 ) = dble( ap( 1 ) )
336  ELSE
337  IF( vl.LT.dble( ap( 1 ) ) .AND. vu.GE.dble( ap( 1 ) ) ) THEN
338  m = 1
339  w( 1 ) = dble( ap( 1 ) )
340  END IF
341  END IF
342  IF( wantz )
343  $ z( 1, 1 ) = cone
344  RETURN
345  END IF
346 *
347 * Get machine constants.
348 *
349  safmin = dlamch( 'Safe minimum' )
350  eps = dlamch( 'Precision' )
351  smlnum = safmin / eps
352  bignum = one / smlnum
353  rmin = sqrt( smlnum )
354  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
355 *
356 * Scale matrix to allowable range, if necessary.
357 *
358  iscale = 0
359  abstll = abstol
360  IF( valeig ) THEN
361  vll = vl
362  vuu = vu
363  ELSE
364  vll = zero
365  vuu = zero
366  END IF
367  anrm = zlanhp( 'M', uplo, n, ap, rwork )
368  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
369  iscale = 1
370  sigma = rmin / anrm
371  ELSE IF( anrm.GT.rmax ) THEN
372  iscale = 1
373  sigma = rmax / anrm
374  END IF
375  IF( iscale.EQ.1 ) THEN
376  CALL zdscal( ( n*( n+1 ) ) / 2, sigma, ap, 1 )
377  IF( abstol.GT.0 )
378  $ abstll = abstol*sigma
379  IF( valeig ) THEN
380  vll = vl*sigma
381  vuu = vu*sigma
382  END IF
383  END IF
384 *
385 * Call ZHPTRD to reduce Hermitian packed matrix to tridiagonal form.
386 *
387  indd = 1
388  inde = indd + n
389  indrwk = inde + n
390  indtau = 1
391  indwrk = indtau + n
392  CALL zhptrd( uplo, n, ap, rwork( indd ), rwork( inde ),
393  $ work( indtau ), iinfo )
394 *
395 * If all eigenvalues are desired and ABSTOL is less than or equal
396 * to zero, then call DSTERF or ZUPGTR and ZSTEQR. If this fails
397 * for some eigenvalue, then try DSTEBZ.
398 *
399  test = .false.
400  IF (indeig) THEN
401  IF (il.EQ.1 .AND. iu.EQ.n) THEN
402  test = .true.
403  END IF
404  END IF
405  IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
406  CALL dcopy( n, rwork( indd ), 1, w, 1 )
407  indee = indrwk + 2*n
408  IF( .NOT.wantz ) THEN
409  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
410  CALL dsterf( n, w, rwork( indee ), info )
411  ELSE
412  CALL zupgtr( uplo, n, ap, work( indtau ), z, ldz,
413  $ work( indwrk ), iinfo )
414  CALL dcopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
415  CALL zsteqr( jobz, n, w, rwork( indee ), z, ldz,
416  $ rwork( indrwk ), info )
417  IF( info.EQ.0 ) THEN
418  DO 10 i = 1, n
419  ifail( i ) = 0
420  10 CONTINUE
421  END IF
422  END IF
423  IF( info.EQ.0 ) THEN
424  m = n
425  GO TO 20
426  END IF
427  info = 0
428  END IF
429 *
430 * Otherwise, call DSTEBZ and, if eigenvectors are desired, ZSTEIN.
431 *
432  IF( wantz ) THEN
433  order = 'B'
434  ELSE
435  order = 'E'
436  END IF
437  indibl = 1
438  indisp = indibl + n
439  indiwk = indisp + n
440  CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
441  $ rwork( indd ), rwork( inde ), m, nsplit, w,
442  $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
443  $ iwork( indiwk ), info )
444 *
445  IF( wantz ) THEN
446  CALL zstein( n, rwork( indd ), rwork( inde ), m, w,
447  $ iwork( indibl ), iwork( indisp ), z, ldz,
448  $ rwork( indrwk ), iwork( indiwk ), ifail, info )
449 *
450 * Apply unitary matrix used in reduction to tridiagonal
451 * form to eigenvectors returned by ZSTEIN.
452 *
453  indwrk = indtau + n
454  CALL zupmtr( 'L', uplo, 'N', n, m, ap, work( indtau ), z, ldz,
455  $ work( indwrk ), iinfo )
456  END IF
457 *
458 * If matrix was scaled, then rescale eigenvalues appropriately.
459 *
460  20 CONTINUE
461  IF( iscale.EQ.1 ) THEN
462  IF( info.EQ.0 ) THEN
463  imax = m
464  ELSE
465  imax = info - 1
466  END IF
467  CALL dscal( imax, one / sigma, w, 1 )
468  END IF
469 *
470 * If eigenvalues are not in order, then sort them, along with
471 * eigenvectors.
472 *
473  IF( wantz ) THEN
474  DO 40 j = 1, m - 1
475  i = 0
476  tmp1 = w( j )
477  DO 30 jj = j + 1, m
478  IF( w( jj ).LT.tmp1 ) THEN
479  i = jj
480  tmp1 = w( jj )
481  END IF
482  30 CONTINUE
483 *
484  IF( i.NE.0 ) THEN
485  itmp1 = iwork( indibl+i-1 )
486  w( i ) = w( j )
487  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
488  w( j ) = tmp1
489  iwork( indibl+j-1 ) = itmp1
490  CALL zswap( n, z( 1, i ), 1, z( 1, j ), 1 )
491  IF( info.NE.0 ) THEN
492  itmp1 = ifail( i )
493  ifail( i ) = ifail( j )
494  ifail( j ) = itmp1
495  END IF
496  END IF
497  40 CONTINUE
498  END IF
499 *
500  RETURN
501 *
502 * End of ZHPEVX
503 *
504  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dstebz(RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, INFO)
DSTEBZ
Definition: dstebz.f:273
subroutine dsterf(N, D, E, INFO)
DSTERF
Definition: dsterf.f:86
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:78
subroutine zsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
ZSTEQR
Definition: zsteqr.f:132
subroutine zupmtr(SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK, INFO)
ZUPMTR
Definition: zupmtr.f:150
subroutine zhptrd(UPLO, N, AP, D, E, TAU, INFO)
ZHPTRD
Definition: zhptrd.f:151
subroutine zstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
ZSTEIN
Definition: zstein.f:182
subroutine zupgtr(UPLO, N, AP, TAU, Q, LDQ, WORK, INFO)
ZUPGTR
Definition: zupgtr.f:114
subroutine zhpevx(JOBZ, RANGE, UPLO, N, AP, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, RWORK, IWORK, IFAIL, INFO)
ZHPEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition: zhpevx.f:240
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79