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

◆ cheevx_2stage()

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

CHEEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for HE matrices

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

Purpose:
 CHEEVX_2STAGE computes selected eigenvalues and, optionally, eigenvectors
 of a complex Hermitian matrix A using the 2stage technique for
 the reduction to tridiagonal.  Eigenvalues and eigenvectors 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.
                  Not available in this release.
[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]A
          A is COMPLEX array, dimension (LDA, N)
          On entry, the Hermitian matrix A.  If UPLO = 'U', the
          leading N-by-N upper triangular part of A contains the
          upper triangular part of the matrix A.  If UPLO = 'L',
          the leading N-by-N lower triangular part of A contains
          the lower triangular part of the matrix A.
          On exit, the lower triangle (if UPLO='L') or the upper
          triangle (if UPLO='U') of A, including the diagonal, is
          destroyed.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[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 A 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)
          On normal exit, the first M elements contain 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 (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The length of the array WORK. LWORK >= 1, when N <= 1;
          otherwise
          If JOBZ = 'N' and N > 1, LWORK must be queried.
                                   LWORK = MAX(1, 8*N, dimension) where
                                   dimension = max(stage1,stage2) + (KD+1)*N + N
                                             = N*KD + N*max(KD+1,FACTOPTNB)
                                               + max(2*KD*KD, KD*NTHREADS)
                                               + (KD+1)*N + N
                                   where KD is the blocking size of the reduction,
                                   FACTOPTNB is the blocking used by the QR or LQ
                                   algorithm, usually FACTOPTNB=128 is a good choice
                                   NTHREADS is the number of threads used when
                                   openMP compilation is enabled, otherwise =1.
          If JOBZ = 'V' and N > 1, LWORK must be queried. Not yet available

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[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.
Further Details:
  All details about the 2stage techniques are available in:

  Azzam Haidar, Hatem Ltaief, and Jack Dongarra.
  Parallel reduction to condensed forms for symmetric eigenvalue problems
  using aggregated fine-grained and memory-aware kernels. In Proceedings
  of 2011 International Conference for High Performance Computing,
  Networking, Storage and Analysis (SC '11), New York, NY, USA,
  Article 8 , 11 pages.
  http://doi.acm.org/10.1145/2063384.2063394

  A. Haidar, J. Kurzak, P. Luszczek, 2013.
  An improved parallel singular value algorithm and its implementation
  for multicore hardware, In Proceedings of 2013 International Conference
  for High Performance Computing, Networking, Storage and Analysis (SC '13).
  Denver, Colorado, USA, 2013.
  Article 90, 12 pages.
  http://doi.acm.org/10.1145/2503210.2503292

  A. Haidar, R. Solca, S. Tomov, T. Schulthess and J. Dongarra.
  A novel hybrid CPU-GPU generalized eigensolver for electronic structure
  calculations based on fine-grained memory aware tasks.
  International Journal of High Performance Computing Applications.
  Volume 28 Issue 2, Pages 196-209, May 2014.
  http://hpc.sagepub.com/content/28/2/196

Definition at line 303 of file cheevx_2stage.f.

306*
307 IMPLICIT NONE
308*
309* -- LAPACK driver routine --
310* -- LAPACK is a software package provided by Univ. of Tennessee, --
311* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
312*
313* .. Scalar Arguments ..
314 CHARACTER JOBZ, RANGE, UPLO
315 INTEGER IL, INFO, IU, LDA, LDZ, LWORK, M, N
316 REAL ABSTOL, VL, VU
317* ..
318* .. Array Arguments ..
319 INTEGER IFAIL( * ), IWORK( * )
320 REAL RWORK( * ), W( * )
321 COMPLEX A( LDA, * ), WORK( * ), Z( LDZ, * )
322* ..
323*
324* =====================================================================
325*
326* .. Parameters ..
327 REAL ZERO, ONE
328 parameter( zero = 0.0e+0, one = 1.0e+0 )
329 COMPLEX CONE
330 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
331* ..
332* .. Local Scalars ..
333 LOGICAL ALLEIG, INDEIG, LOWER, LQUERY, TEST, VALEIG,
334 $ WANTZ
335 CHARACTER ORDER
336 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
337 $ INDISP, INDIWK, INDRWK, INDTAU, INDWRK, ISCALE,
338 $ ITMP1, J, JJ, LLWORK,
339 $ NSPLIT, LWMIN, LHTRD, LWTRD, KD, IB, INDHOUS
340 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
341 $ SIGMA, SMLNUM, TMP1, VLL, VUU
342* ..
343* .. External Functions ..
344 LOGICAL LSAME
345 INTEGER ILAENV2STAGE
346 REAL SLAMCH, CLANHE, SROUNDUP_LWORK
347 EXTERNAL lsame, slamch, clanhe, ilaenv2stage,
349* ..
350* .. External Subroutines ..
351 EXTERNAL scopy, sscal, sstebz, ssterf, xerbla, csscal,
354* ..
355* .. Intrinsic Functions ..
356 INTRINSIC real, max, min, sqrt
357* ..
358* .. Executable Statements ..
359*
360* Test the input parameters.
361*
362 lower = lsame( uplo, 'L' )
363 wantz = lsame( jobz, 'V' )
364 alleig = lsame( range, 'A' )
365 valeig = lsame( range, 'V' )
366 indeig = lsame( range, 'I' )
367 lquery = ( lwork.EQ.-1 )
368*
369 info = 0
370 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
371 info = -1
372 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
373 info = -2
374 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
375 info = -3
376 ELSE IF( n.LT.0 ) THEN
377 info = -4
378 ELSE IF( lda.LT.max( 1, n ) ) THEN
379 info = -6
380 ELSE
381 IF( valeig ) THEN
382 IF( n.GT.0 .AND. vu.LE.vl )
383 $ info = -8
384 ELSE IF( indeig ) THEN
385 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
386 info = -9
387 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
388 info = -10
389 END IF
390 END IF
391 END IF
392 IF( info.EQ.0 ) THEN
393 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) ) THEN
394 info = -15
395 END IF
396 END IF
397*
398 IF( info.EQ.0 ) THEN
399 IF( n.LE.1 ) THEN
400 lwmin = 1
401 work( 1 ) = sroundup_lwork(lwmin)
402 ELSE
403 kd = ilaenv2stage( 1, 'CHETRD_2STAGE', jobz,
404 $ n, -1, -1, -1 )
405 ib = ilaenv2stage( 2, 'CHETRD_2STAGE', jobz,
406 $ n, kd, -1, -1 )
407 lhtrd = ilaenv2stage( 3, 'CHETRD_2STAGE', jobz,
408 $ n, kd, ib, -1 )
409 lwtrd = ilaenv2stage( 4, 'CHETRD_2STAGE', jobz,
410 $ n, kd, ib, -1 )
411 lwmin = n + lhtrd + lwtrd
412 work( 1 ) = sroundup_lwork(lwmin)
413 END IF
414*
415 IF( lwork.LT.lwmin .AND. .NOT.lquery )
416 $ info = -17
417 END IF
418*
419 IF( info.NE.0 ) THEN
420 CALL xerbla( 'CHEEVX_2STAGE', -info )
421 RETURN
422 ELSE IF( lquery ) THEN
423 RETURN
424 END IF
425*
426* Quick return if possible
427*
428 m = 0
429 IF( n.EQ.0 ) THEN
430 RETURN
431 END IF
432*
433 IF( n.EQ.1 ) THEN
434 IF( alleig .OR. indeig ) THEN
435 m = 1
436 w( 1 ) = real( a( 1, 1 ) )
437 ELSE IF( valeig ) THEN
438 IF( vl.LT.real( a( 1, 1 ) ) .AND. vu.GE.real( a( 1, 1 ) ) )
439 $ THEN
440 m = 1
441 w( 1 ) = real( a( 1, 1 ) )
442 END IF
443 END IF
444 IF( wantz )
445 $ z( 1, 1 ) = cone
446 RETURN
447 END IF
448*
449* Get machine constants.
450*
451 safmin = slamch( 'Safe minimum' )
452 eps = slamch( 'Precision' )
453 smlnum = safmin / eps
454 bignum = one / smlnum
455 rmin = sqrt( smlnum )
456 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
457*
458* Scale matrix to allowable range, if necessary.
459*
460 iscale = 0
461 abstll = abstol
462 IF( valeig ) THEN
463 vll = vl
464 vuu = vu
465 END IF
466 anrm = clanhe( 'M', uplo, n, a, lda, rwork )
467 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
468 iscale = 1
469 sigma = rmin / anrm
470 ELSE IF( anrm.GT.rmax ) THEN
471 iscale = 1
472 sigma = rmax / anrm
473 END IF
474 IF( iscale.EQ.1 ) THEN
475 IF( lower ) THEN
476 DO 10 j = 1, n
477 CALL csscal( n-j+1, sigma, a( j, j ), 1 )
478 10 CONTINUE
479 ELSE
480 DO 20 j = 1, n
481 CALL csscal( j, sigma, a( 1, j ), 1 )
482 20 CONTINUE
483 END IF
484 IF( abstol.GT.0 )
485 $ abstll = abstol*sigma
486 IF( valeig ) THEN
487 vll = vl*sigma
488 vuu = vu*sigma
489 END IF
490 END IF
491*
492* Call CHETRD_2STAGE to reduce Hermitian matrix to tridiagonal form.
493*
494 indd = 1
495 inde = indd + n
496 indrwk = inde + n
497 indtau = 1
498 indhous = indtau + n
499 indwrk = indhous + lhtrd
500 llwork = lwork - indwrk + 1
501*
502 CALL chetrd_2stage( jobz, uplo, n, a, lda, rwork( indd ),
503 $ rwork( inde ), work( indtau ),
504 $ work( indhous ), lhtrd, work( indwrk ),
505 $ llwork, iinfo )
506*
507* If all eigenvalues are desired and ABSTOL is less than or equal to
508* zero, then call SSTERF or CUNGTR and CSTEQR. If this fails for
509* some eigenvalue, then try SSTEBZ.
510*
511 test = .false.
512 IF( indeig ) THEN
513 IF( il.EQ.1 .AND. iu.EQ.n ) THEN
514 test = .true.
515 END IF
516 END IF
517 IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
518 CALL scopy( n, rwork( indd ), 1, w, 1 )
519 indee = indrwk + 2*n
520 IF( .NOT.wantz ) THEN
521 CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
522 CALL ssterf( n, w, rwork( indee ), info )
523 ELSE
524 CALL clacpy( 'A', n, n, a, lda, z, ldz )
525 CALL cungtr( uplo, n, z, ldz, work( indtau ),
526 $ work( indwrk ), llwork, iinfo )
527 CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
528 CALL csteqr( jobz, n, w, rwork( indee ), z, ldz,
529 $ rwork( indrwk ), info )
530 IF( info.EQ.0 ) THEN
531 DO 30 i = 1, n
532 ifail( i ) = 0
533 30 CONTINUE
534 END IF
535 END IF
536 IF( info.EQ.0 ) THEN
537 m = n
538 GO TO 40
539 END IF
540 info = 0
541 END IF
542*
543* Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN.
544*
545 IF( wantz ) THEN
546 order = 'B'
547 ELSE
548 order = 'E'
549 END IF
550 indibl = 1
551 indisp = indibl + n
552 indiwk = indisp + n
553 CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
554 $ rwork( indd ), rwork( inde ), m, nsplit, w,
555 $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
556 $ iwork( indiwk ), info )
557*
558 IF( wantz ) THEN
559 CALL cstein( n, rwork( indd ), rwork( inde ), m, w,
560 $ iwork( indibl ), iwork( indisp ), z, ldz,
561 $ rwork( indrwk ), iwork( indiwk ), ifail, info )
562*
563* Apply unitary matrix used in reduction to tridiagonal
564* form to eigenvectors returned by CSTEIN.
565*
566 CALL cunmtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
567 $ ldz, work( indwrk ), llwork, iinfo )
568 END IF
569*
570* If matrix was scaled, then rescale eigenvalues appropriately.
571*
572 40 CONTINUE
573 IF( iscale.EQ.1 ) THEN
574 IF( info.EQ.0 ) THEN
575 imax = m
576 ELSE
577 imax = info - 1
578 END IF
579 CALL sscal( imax, one / sigma, w, 1 )
580 END IF
581*
582* If eigenvalues are not in order, then sort them, along with
583* eigenvectors.
584*
585 IF( wantz ) THEN
586 DO 60 j = 1, m - 1
587 i = 0
588 tmp1 = w( j )
589 DO 50 jj = j + 1, m
590 IF( w( jj ).LT.tmp1 ) THEN
591 i = jj
592 tmp1 = w( jj )
593 END IF
594 50 CONTINUE
595*
596 IF( i.NE.0 ) THEN
597 itmp1 = iwork( indibl+i-1 )
598 w( i ) = w( j )
599 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
600 w( j ) = tmp1
601 iwork( indibl+j-1 ) = itmp1
602 CALL cswap( n, z( 1, i ), 1, z( 1, j ), 1 )
603 IF( info.NE.0 ) THEN
604 itmp1 = ifail( i )
605 ifail( i ) = ifail( j )
606 ifail( j ) = itmp1
607 END IF
608 END IF
609 60 CONTINUE
610 END IF
611*
612* Set WORK(1) to optimal complex workspace size.
613*
614 work( 1 ) = sroundup_lwork(lwmin)
615*
616 RETURN
617*
618* End of CHEEVX_2STAGE
619*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine chetrd_2stage(vect, uplo, n, a, lda, d, e, tau, hous2, lhous2, work, lwork, info)
CHETRD_2STAGE
integer function ilaenv2stage(ispec, name, opts, n1, n2, n3, n4)
ILAENV2STAGE
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
real function clanhe(norm, uplo, n, a, lda, work)
CLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition clanhe.f:124
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
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 cungtr(uplo, n, a, lda, tau, work, lwork, info)
CUNGTR
Definition cungtr.f:123
subroutine cunmtr(side, uplo, trans, m, n, a, lda, tau, c, ldc, work, lwork, info)
CUNMTR
Definition cunmtr.f:172
Here is the call graph for this function:
Here is the caller graph for this function: