LAPACK 3.12.0
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 "Computing Small Singular Values of Bidiagonal Matrices
          with Guaranteed High Relative Accuracy," 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 237 of file chpevx.f.

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 REAL ABSTOL, VL, VU
249* ..
250* .. Array Arguments ..
251 INTEGER IFAIL( * ), IWORK( * )
252 REAL RWORK( * ), W( * )
253 COMPLEX AP( * ), WORK( * ), Z( LDZ, * )
254* ..
255*
256* =====================================================================
257*
258* .. Parameters ..
259 REAL ZERO, ONE
260 parameter( zero = 0.0e0, one = 1.0e0 )
261 COMPLEX CONE
262 parameter( cone = ( 1.0e0, 0.0e0 ) )
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 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
271 $ SIGMA, SMLNUM, TMP1, VLL, VUU
272* ..
273* .. External Functions ..
274 LOGICAL LSAME
275 REAL CLANHP, SLAMCH
276 EXTERNAL lsame, clanhp, slamch
277* ..
278* .. External Subroutines ..
279 EXTERNAL chptrd, csscal, cstein, csteqr, cswap, cupgtr,
281* ..
282* .. Intrinsic Functions ..
283 INTRINSIC max, min, real, 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( '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, 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 sscal( 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 cswap( 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 CHPEVX
502*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine chptrd(uplo, n, ap, d, e, tau, info)
CHPTRD
Definition chptrd.f:151
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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:117
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine sstebz(range, order, n, vl, vu, il, iu, abstol, d, e, m, nsplit, w, iblock, isplit, work, iwork, info)
SSTEBZ
Definition sstebz.f:273
subroutine cstein(n, d, e, m, w, iblock, isplit, z, ldz, work, iwork, ifail, info)
CSTEIN
Definition cstein.f:182
subroutine csteqr(compz, n, d, e, z, ldz, work, info)
CSTEQR
Definition csteqr.f:132
subroutine ssterf(n, d, e, info)
SSTERF
Definition ssterf.f:86
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
subroutine cupgtr(uplo, n, ap, tau, q, ldq, work, info)
CUPGTR
Definition cupgtr.f:114
subroutine cupmtr(side, uplo, trans, m, n, ap, tau, c, ldc, work, info)
CUPMTR
Definition cupmtr.f:150
Here is the call graph for this function:
Here is the caller graph for this function: