LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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 hpevx
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,
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 indisp = 1 + n
438 indiwk = indisp + n
439 CALL dstebz( range, order, n, vll, vuu, il, iu, abstll,
440 $ rwork( indd ), rwork( inde ), m, nsplit, w,
441 $ iwork( 1 ), iwork( indisp ), rwork( indrwk ),
442 $ iwork( indiwk ), info )
443*
444 IF( wantz ) THEN
445 CALL zstein( n, rwork( indd ), rwork( inde ), m, w,
446 $ iwork( 1 ), iwork( indisp ), z, ldz,
447 $ rwork( indrwk ), iwork( indiwk ), ifail, info )
448*
449* Apply unitary matrix used in reduction to tridiagonal
450* form to eigenvectors returned by ZSTEIN.
451*
452 indwrk = indtau + n
453 CALL zupmtr( 'L', uplo, 'N', n, m, ap, work( indtau ), z, ldz,
454 $ work( indwrk ), iinfo )
455 END IF
456*
457* If matrix was scaled, then rescale eigenvalues appropriately.
458*
459 20 CONTINUE
460 IF( iscale.EQ.1 ) THEN
461 IF( info.EQ.0 ) THEN
462 imax = m
463 ELSE
464 imax = info - 1
465 END IF
466 CALL dscal( imax, one / sigma, w, 1 )
467 END IF
468*
469* If eigenvalues are not in order, then sort them, along with
470* eigenvectors.
471*
472 IF( wantz ) THEN
473 DO 40 j = 1, m - 1
474 i = 0
475 tmp1 = w( j )
476 DO 30 jj = j + 1, m
477 IF( w( jj ).LT.tmp1 ) THEN
478 i = jj
479 tmp1 = w( jj )
480 END IF
481 30 CONTINUE
482*
483 IF( i.NE.0 ) THEN
484 itmp1 = iwork( 1 + i-1 )
485 w( i ) = w( j )
486 iwork( 1 + i-1 ) = iwork( 1 + j-1 )
487 w( j ) = tmp1
488 iwork( 1 + j-1 ) = itmp1
489 CALL zswap( n, z( 1, i ), 1, z( 1, j ), 1 )
490 IF( info.NE.0 ) THEN
491 itmp1 = ifail( i )
492 ifail( i ) = ifail( j )
493 ifail( j ) = itmp1
494 END IF
495 END IF
496 40 CONTINUE
497 END IF
498*
499 RETURN
500*
501* End of ZHPEVX
502*
503 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
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 zhptrd(uplo, n, ap, d, e, tau, info)
ZHPTRD
Definition zhptrd.f:151
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
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 zstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
ZSTEIN
Definition zstein.f:182
subroutine zsteqr(compz, n, d, e, z, ldz, work, info)
ZSTEQR
Definition zsteqr.f:132
subroutine dsterf(n, d, e, info)
DSTERF
Definition dsterf.f:86
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
subroutine zupgtr(uplo, n, ap, tau, q, ldq, work, info)
ZUPGTR
Definition zupgtr.f:114
subroutine zupmtr(side, uplo, trans, m, n, ap, tau, c, ldc, work, info)
ZUPMTR
Definition zupmtr.f:150