LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ chpevx()

subroutine chpevx ( character jobz,
character range,
character uplo,
integer n,
complex, dimension( * ) ap,
real vl,
real vu,
integer il,
integer iu,
real abstol,
integer m,
real, dimension( * ) w,
complex, dimension( ldz, * ) z,
integer ldz,
complex, dimension( * ) work,
real, dimension( * ) rwork,
integer, dimension( * ) iwork,
integer, dimension( * ) ifail,
integer info )

CHPEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrices

Download CHPEVX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> CHPEVX computes selected eigenvalues and, optionally, eigenvectors
!> of a complex Hermitian matrix A in packed storage.
!> Eigenvalues/vectors can be selected by specifying either a range of
!> values or a range of indices for the desired eigenvalues.
!>
Parameters
[in]JOBZ
!>          JOBZ is CHARACTER*1
!>          = 'N':  Compute eigenvalues only;
!>          = 'V':  Compute eigenvalues and eigenvectors.
!>
[in]RANGE
!>          RANGE is CHARACTER*1
!>          = 'A': all eigenvalues will be found;
!>          = 'V': all eigenvalues in the half-open interval (VL,VU]
!>                 will be found;
!>          = 'I': the IL-th through IU-th eigenvalues will be found.
!>
[in]UPLO
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!>
[in]N
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!>
[in,out]AP
!>          AP is COMPLEX array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the Hermitian matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>
!>          On exit, AP is overwritten by values generated during the
!>          reduction to tridiagonal form.  If UPLO = 'U', the diagonal
!>          and first superdiagonal of the tridiagonal matrix T overwrite
!>          the corresponding elements of A, and if UPLO = 'L', the
!>          diagonal and first subdiagonal of T overwrite the
!>          corresponding elements of A.
!>
[in]VL
!>          VL is REAL
!>          If RANGE='V', the lower bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]VU
!>          VU is REAL
!>          If RANGE='V', the upper bound of the interval to
!>          be searched for eigenvalues. VL < VU.
!>          Not referenced if RANGE = 'A' or 'I'.
!>
[in]IL
!>          IL is INTEGER
!>          If RANGE='I', the index of the
!>          smallest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]IU
!>          IU is INTEGER
!>          If RANGE='I', the index of the
!>          largest eigenvalue to be returned.
!>          1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0.
!>          Not referenced if RANGE = 'A' or 'V'.
!>
[in]ABSTOL
!>          ABSTOL is REAL
!>          The absolute error tolerance for the eigenvalues.
!>          An approximate eigenvalue is accepted as converged
!>          when it is determined to lie in an interval [a,b]
!>          of width less than or equal to
!>
!>                  ABSTOL + EPS *   max( |a|,|b| ) ,
!>
!>          where EPS is the machine precision.  If ABSTOL is less than
!>          or equal to zero, then  EPS*|T|  will be used in its place,
!>          where |T| is the 1-norm of the tridiagonal matrix obtained
!>          by reducing AP to tridiagonal form.
!>
!>          Eigenvalues will be computed most accurately when ABSTOL is
!>          set to twice the underflow threshold 2*SLAMCH('S'), not zero.
!>          If this routine returns with INFO>0, indicating that some
!>          eigenvectors did not converge, try setting ABSTOL to
!>          2*SLAMCH('S').
!>
!>          See  by Demmel and
!>          Kahan, LAPACK Working Note #3.
!>
[out]M
!>          M is INTEGER
!>          The total number of eigenvalues found.  0 <= M <= N.
!>          If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1.
!>
[out]W
!>          W is REAL array, dimension (N)
!>          If INFO = 0, the selected eigenvalues in ascending order.
!>
[out]Z
!>          Z is COMPLEX array, dimension (LDZ, max(1,M))
!>          If JOBZ = 'V', then if INFO = 0, the first M columns of Z
!>          contain the orthonormal eigenvectors of the matrix A
!>          corresponding to the selected eigenvalues, with the i-th
!>          column of Z holding the eigenvector associated with W(i).
!>          If an eigenvector fails to converge, then that column of Z
!>          contains the latest approximation to the eigenvector, and
!>          the index of the eigenvector is returned in IFAIL.
!>          If JOBZ = 'N', then Z is not referenced.
!>          Note: the user must ensure that at least max(1,M) columns are
!>          supplied in the array Z; if RANGE = 'V', the exact value of M
!>          is not known in advance and an upper bound must be used.
!>
[in]LDZ
!>          LDZ is INTEGER
!>          The leading dimension of the array Z.  LDZ >= 1, and if
!>          JOBZ = 'V', LDZ >= max(1,N).
!>
[out]WORK
!>          WORK is COMPLEX array, dimension (2*N)
!>
[out]RWORK
!>          RWORK is REAL array, dimension (7*N)
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (5*N)
!>
[out]IFAIL
!>          IFAIL is INTEGER array, dimension (N)
!>          If JOBZ = 'V', then if INFO = 0, the first M elements of
!>          IFAIL are zero.  If INFO > 0, then IFAIL contains the
!>          indices of the eigenvectors that failed to converge.
!>          If JOBZ = 'N', then IFAIL is not referenced.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!>          > 0:  if INFO = i, then i eigenvectors failed to converge.
!>                Their indices are stored in array IFAIL.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 235 of file chpevx.f.

238*
239* -- LAPACK driver routine --
240* -- LAPACK is a software package provided by Univ. of Tennessee, --
241* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
242*
243* .. Scalar Arguments ..
244 CHARACTER JOBZ, RANGE, UPLO
245 INTEGER IL, INFO, IU, LDZ, M, N
246 REAL ABSTOL, VL, VU
247* ..
248* .. Array Arguments ..
249 INTEGER IFAIL( * ), IWORK( * )
250 REAL RWORK( * ), W( * )
251 COMPLEX AP( * ), WORK( * ), Z( LDZ, * )
252* ..
253*
254* =====================================================================
255*
256* .. Parameters ..
257 REAL ZERO, ONE
258 parameter( zero = 0.0e0, one = 1.0e0 )
259 COMPLEX CONE
260 parameter( cone = ( 1.0e0, 0.0e0 ) )
261* ..
262* .. Local Scalars ..
263 LOGICAL ALLEIG, INDEIG, TEST, VALEIG, WANTZ
264 CHARACTER ORDER
265 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE,
266 $ INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE,
267 $ ITMP1, J, JJ, NSPLIT
268 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
269 $ SIGMA, SMLNUM, TMP1, VLL, VUU
270* ..
271* .. External Functions ..
272 LOGICAL LSAME
273 REAL CLANHP, SLAMCH
274 EXTERNAL lsame, clanhp, slamch
275* ..
276* .. External Subroutines ..
277 EXTERNAL chptrd, csscal, cstein, csteqr, cswap,
278 $ cupgtr,
279 $ cupmtr, scopy, sscal, sstebz, ssterf, xerbla
280* ..
281* .. Intrinsic Functions ..
282 INTRINSIC max, min, real, sqrt
283* ..
284* .. Executable Statements ..
285*
286* Test the input parameters.
287*
288 wantz = lsame( jobz, 'V' )
289 alleig = lsame( range, 'A' )
290 valeig = lsame( range, 'V' )
291 indeig = lsame( range, 'I' )
292*
293 info = 0
294 IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
295 info = -1
296 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
297 info = -2
298 ELSE IF( .NOT.( lsame( uplo, 'L' ) .OR.
299 $ 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( 'CHPEVX', -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 ) = real( ap( 1 ) )
336 ELSE
337 IF( vl.LT.real( ap( 1 ) ) .AND. vu.GE.real( ap( 1 ) ) ) THEN
338 m = 1
339 w( 1 ) = real( 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 = slamch( 'Safe minimum' )
350 eps = slamch( '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 ENDIF
367 anrm = clanhp( '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 csscal( ( 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 CHPTRD 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 chptrd( 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 SSTERF or CUPGTR and CSTEQR. If this fails
397* for some eigenvalue, then try SSTEBZ.
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 scopy( n, rwork( indd ), 1, w, 1 )
407 indee = indrwk + 2*n
408 IF( .NOT.wantz ) THEN
409 CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
410 CALL ssterf( n, w, rwork( indee ), info )
411 ELSE
412 CALL cupgtr( uplo, n, ap, work( indtau ), z, ldz,
413 $ work( indwrk ), iinfo )
414 CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
415 CALL csteqr( 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 SSTEBZ and, if eigenvectors are desired, CSTEIN.
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 sstebz( 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 cstein( 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 CSTEIN.
451*
452 indwrk = indtau + n
453 CALL cupmtr( 'L', uplo, 'N', n, m, ap, work( indtau ), z,
454 $ 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 sscal( 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( 1 + i-1 )
486 w( i ) = w( j )
487 iwork( 1 + i-1 ) = iwork( 1 + j-1 )
488 w( j ) = tmp1
489 iwork( 1 + j-1 ) = itmp1
490 CALL cswap( 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 CHPEVX
503*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
chptrd
subroutine chptrd(uplo, n, ap, d, e, tau, info)
CHPTRD
Definition chptrd.f:149
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
clanhp
real function clanhp(norm, uplo, n, ap, work)
CLANHP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition clanhp.f:115
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
sstebz
subroutine sstebz(range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, iwork, info)
SSTEBZ
Definition sstebz.f:272
cstein
subroutine cstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
CSTEIN
Definition cstein.f:180
csteqr
subroutine csteqr(compz, n, d, e, z, ldz, work, info)
CSTEQR
Definition csteqr.f:130
ssterf
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:84
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
cupgtr
subroutine cupgtr(uplo, n, ap, tau, q, ldq, work, info)
CUPGTR
Definition cupgtr.f:112
cupmtr
subroutine cupmtr(side, uplo, trans, m, n, ap, tau, c, ldc, work, info)
CUPMTR
Definition cupmtr.f:149
Here is the call graph for this function:
Here is the caller graph for this function: