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

◆ chbevx_2stage()

subroutine chbevx_2stage ( character  jobz,
character  range,
character  uplo,
integer  n,
integer  kd,
complex, dimension( ldab, * )  ab,
integer  ldab,
complex, dimension( ldq, * )  q,
integer  ldq,
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 
)

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

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

Purpose:
 CHBEVX_2STAGE computes selected eigenvalues and, optionally, eigenvectors
 of a complex Hermitian band 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]KD
          KD is INTEGER
          The number of superdiagonals of the matrix A if UPLO = 'U',
          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
[in,out]AB
          AB is COMPLEX array, dimension (LDAB, N)
          On entry, the upper or lower triangle of the Hermitian band
          matrix A, stored in the first KD+1 rows of the array.  The
          j-th column of A is stored in the j-th column of the array AB
          as follows:
          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).

          On exit, AB is overwritten by values generated during the
          reduction to tridiagonal form.
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array AB.  LDAB >= KD + 1.
[out]Q
          Q is COMPLEX array, dimension (LDQ, N)
          If JOBZ = 'V', the N-by-N unitary matrix used in the
                          reduction to tridiagonal form.
          If JOBZ = 'N', the array Q is not referenced.
[in]LDQ
          LDQ is INTEGER
          The leading dimension of the array Q.  If JOBZ = 'V', then
          LDQ >= 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 AB 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)
          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 (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, dimension) where
                                   dimension = (2KD+1)*N + KD*NTHREADS
                                   where KD is the size of the band.
                                   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 sizes of the WORK, RWORK and
          IWORK arrays, returns these values as the first entries of
          the WORK, RWORK and IWORK arrays, and no error message
          related to LWORK or LRWORK or LIWORK 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 323 of file chbevx_2stage.f.

327*
328 IMPLICIT NONE
329*
330* -- LAPACK driver routine --
331* -- LAPACK is a software package provided by Univ. of Tennessee, --
332* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
333*
334* .. Scalar Arguments ..
335 CHARACTER JOBZ, RANGE, UPLO
336 INTEGER IL, INFO, IU, KD, LDAB, LDQ, LDZ, M, N, LWORK
337 REAL ABSTOL, VL, VU
338* ..
339* .. Array Arguments ..
340 INTEGER IFAIL( * ), IWORK( * )
341 REAL RWORK( * ), W( * )
342 COMPLEX AB( LDAB, * ), Q( LDQ, * ), WORK( * ),
343 $ Z( LDZ, * )
344* ..
345*
346* =====================================================================
347*
348* .. Parameters ..
349 REAL ZERO, ONE
350 parameter( zero = 0.0e0, one = 1.0e0 )
351 COMPLEX CZERO, CONE
352 parameter( czero = ( 0.0e0, 0.0e0 ),
353 $ cone = ( 1.0e0, 0.0e0 ) )
354* ..
355* .. Local Scalars ..
356 LOGICAL ALLEIG, INDEIG, LOWER, TEST, VALEIG, WANTZ,
357 $ LQUERY
358 CHARACTER ORDER
359 INTEGER I, IINFO, IMAX, INDD, INDE, INDEE, INDIBL,
360 $ INDISP, INDIWK, INDRWK, INDWRK, ISCALE, ITMP1,
361 $ LLWORK, LWMIN, LHTRD, LWTRD, IB, INDHOUS,
362 $ J, JJ, NSPLIT
363 REAL ABSTLL, ANRM, BIGNUM, EPS, RMAX, RMIN, SAFMIN,
364 $ SIGMA, SMLNUM, TMP1, VLL, VUU
365 COMPLEX CTMP1
366* ..
367* .. External Functions ..
368 LOGICAL LSAME
369 INTEGER ILAENV2STAGE
370 REAL SLAMCH, CLANHB, SROUNDUP_LWORK
371 EXTERNAL lsame, slamch, clanhb, ilaenv2stage,
373* ..
374* .. External Subroutines ..
375 EXTERNAL scopy, sscal, sstebz, ssterf, xerbla, ccopy,
378* ..
379* .. Intrinsic Functions ..
380 INTRINSIC real, max, min, sqrt
381* ..
382* .. Executable Statements ..
383*
384* Test the input parameters.
385*
386 wantz = lsame( jobz, 'V' )
387 alleig = lsame( range, 'A' )
388 valeig = lsame( range, 'V' )
389 indeig = lsame( range, 'I' )
390 lower = lsame( uplo, 'L' )
391 lquery = ( lwork.EQ.-1 )
392*
393 info = 0
394 IF( .NOT.( lsame( jobz, 'N' ) ) ) THEN
395 info = -1
396 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) ) THEN
397 info = -2
398 ELSE IF( .NOT.( lower .OR. lsame( uplo, 'U' ) ) ) THEN
399 info = -3
400 ELSE IF( n.LT.0 ) THEN
401 info = -4
402 ELSE IF( kd.LT.0 ) THEN
403 info = -5
404 ELSE IF( ldab.LT.kd+1 ) THEN
405 info = -7
406 ELSE IF( wantz .AND. ldq.LT.max( 1, n ) ) THEN
407 info = -9
408 ELSE
409 IF( valeig ) THEN
410 IF( n.GT.0 .AND. vu.LE.vl )
411 $ info = -11
412 ELSE IF( indeig ) THEN
413 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
414 info = -12
415 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
416 info = -13
417 END IF
418 END IF
419 END IF
420 IF( info.EQ.0 ) THEN
421 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n ) )
422 $ info = -18
423 END IF
424*
425 IF( info.EQ.0 ) THEN
426 IF( n.LE.1 ) THEN
427 lwmin = 1
428 work( 1 ) = sroundup_lwork(lwmin)
429 ELSE
430 ib = ilaenv2stage( 2, 'CHETRD_HB2ST', jobz,
431 $ n, kd, -1, -1 )
432 lhtrd = ilaenv2stage( 3, 'CHETRD_HB2ST', jobz,
433 $ n, kd, ib, -1 )
434 lwtrd = ilaenv2stage( 4, 'CHETRD_HB2ST', jobz,
435 $ n, kd, ib, -1 )
436 lwmin = lhtrd + lwtrd
437 work( 1 ) = sroundup_lwork(lwmin)
438 ENDIF
439*
440 IF( lwork.LT.lwmin .AND. .NOT.lquery )
441 $ info = -20
442 END IF
443*
444 IF( info.NE.0 ) THEN
445 CALL xerbla( 'CHBEVX_2STAGE', -info )
446 RETURN
447 ELSE IF( lquery ) THEN
448 RETURN
449 END IF
450*
451* Quick return if possible
452*
453 m = 0
454 IF( n.EQ.0 )
455 $ RETURN
456*
457 IF( n.EQ.1 ) THEN
458 m = 1
459 IF( lower ) THEN
460 ctmp1 = ab( 1, 1 )
461 ELSE
462 ctmp1 = ab( kd+1, 1 )
463 END IF
464 tmp1 = real( ctmp1 )
465 IF( valeig ) THEN
466 IF( .NOT.( vl.LT.tmp1 .AND. vu.GE.tmp1 ) )
467 $ m = 0
468 END IF
469 IF( m.EQ.1 ) THEN
470 w( 1 ) = real( ctmp1 )
471 IF( wantz )
472 $ z( 1, 1 ) = cone
473 END IF
474 RETURN
475 END IF
476*
477* Get machine constants.
478*
479 safmin = slamch( 'Safe minimum' )
480 eps = slamch( 'Precision' )
481 smlnum = safmin / eps
482 bignum = one / smlnum
483 rmin = sqrt( smlnum )
484 rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
485*
486* Scale matrix to allowable range, if necessary.
487*
488 iscale = 0
489 abstll = abstol
490 IF( valeig ) THEN
491 vll = vl
492 vuu = vu
493 ELSE
494 vll = zero
495 vuu = zero
496 END IF
497 anrm = clanhb( 'M', uplo, n, kd, ab, ldab, rwork )
498 IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
499 iscale = 1
500 sigma = rmin / anrm
501 ELSE IF( anrm.GT.rmax ) THEN
502 iscale = 1
503 sigma = rmax / anrm
504 END IF
505 IF( iscale.EQ.1 ) THEN
506 IF( lower ) THEN
507 CALL clascl( 'B', kd, kd, one, sigma, n, n, ab, ldab, info )
508 ELSE
509 CALL clascl( 'Q', kd, kd, one, sigma, n, n, ab, ldab, info )
510 END IF
511 IF( abstol.GT.0 )
512 $ abstll = abstol*sigma
513 IF( valeig ) THEN
514 vll = vl*sigma
515 vuu = vu*sigma
516 END IF
517 END IF
518*
519* Call CHBTRD_HB2ST to reduce Hermitian band matrix to tridiagonal form.
520*
521 indd = 1
522 inde = indd + n
523 indrwk = inde + n
524*
525 indhous = 1
526 indwrk = indhous + lhtrd
527 llwork = lwork - indwrk + 1
528*
529 CALL chetrd_hb2st( 'N', jobz, uplo, n, kd, ab, ldab,
530 $ rwork( indd ), rwork( inde ), work( indhous ),
531 $ lhtrd, work( indwrk ), llwork, iinfo )
532*
533* If all eigenvalues are desired and ABSTOL is less than or equal
534* to zero, then call SSTERF or CSTEQR. If this fails for some
535* eigenvalue, then try SSTEBZ.
536*
537 test = .false.
538 IF (indeig) THEN
539 IF (il.EQ.1 .AND. iu.EQ.n) THEN
540 test = .true.
541 END IF
542 END IF
543 IF ((alleig .OR. test) .AND. (abstol.LE.zero)) THEN
544 CALL scopy( n, rwork( indd ), 1, w, 1 )
545 indee = indrwk + 2*n
546 IF( .NOT.wantz ) THEN
547 CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
548 CALL ssterf( n, w, rwork( indee ), info )
549 ELSE
550 CALL clacpy( 'A', n, n, q, ldq, z, ldz )
551 CALL scopy( n-1, rwork( inde ), 1, rwork( indee ), 1 )
552 CALL csteqr( jobz, n, w, rwork( indee ), z, ldz,
553 $ rwork( indrwk ), info )
554 IF( info.EQ.0 ) THEN
555 DO 10 i = 1, n
556 ifail( i ) = 0
557 10 CONTINUE
558 END IF
559 END IF
560 IF( info.EQ.0 ) THEN
561 m = n
562 GO TO 30
563 END IF
564 info = 0
565 END IF
566*
567* Otherwise, call SSTEBZ and, if eigenvectors are desired, CSTEIN.
568*
569 IF( wantz ) THEN
570 order = 'B'
571 ELSE
572 order = 'E'
573 END IF
574 indibl = 1
575 indisp = indibl + n
576 indiwk = indisp + n
577 CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
578 $ rwork( indd ), rwork( inde ), m, nsplit, w,
579 $ iwork( indibl ), iwork( indisp ), rwork( indrwk ),
580 $ iwork( indiwk ), info )
581*
582 IF( wantz ) THEN
583 CALL cstein( n, rwork( indd ), rwork( inde ), m, w,
584 $ iwork( indibl ), iwork( indisp ), z, ldz,
585 $ rwork( indrwk ), iwork( indiwk ), ifail, info )
586*
587* Apply unitary matrix used in reduction to tridiagonal
588* form to eigenvectors returned by CSTEIN.
589*
590 DO 20 j = 1, m
591 CALL ccopy( n, z( 1, j ), 1, work( 1 ), 1 )
592 CALL cgemv( 'N', n, n, cone, q, ldq, work, 1, czero,
593 $ z( 1, j ), 1 )
594 20 CONTINUE
595 END IF
596*
597* If matrix was scaled, then rescale eigenvalues appropriately.
598*
599 30 CONTINUE
600 IF( iscale.EQ.1 ) THEN
601 IF( info.EQ.0 ) THEN
602 imax = m
603 ELSE
604 imax = info - 1
605 END IF
606 CALL sscal( imax, one / sigma, w, 1 )
607 END IF
608*
609* If eigenvalues are not in order, then sort them, along with
610* eigenvectors.
611*
612 IF( wantz ) THEN
613 DO 50 j = 1, m - 1
614 i = 0
615 tmp1 = w( j )
616 DO 40 jj = j + 1, m
617 IF( w( jj ).LT.tmp1 ) THEN
618 i = jj
619 tmp1 = w( jj )
620 END IF
621 40 CONTINUE
622*
623 IF( i.NE.0 ) THEN
624 itmp1 = iwork( indibl+i-1 )
625 w( i ) = w( j )
626 iwork( indibl+i-1 ) = iwork( indibl+j-1 )
627 w( j ) = tmp1
628 iwork( indibl+j-1 ) = itmp1
629 CALL cswap( n, z( 1, i ), 1, z( 1, j ), 1 )
630 IF( info.NE.0 ) THEN
631 itmp1 = ifail( i )
632 ifail( i ) = ifail( j )
633 ifail( j ) = itmp1
634 END IF
635 END IF
636 50 CONTINUE
637 END IF
638*
639* Set WORK(1) to optimal workspace size.
640*
641 work( 1 ) = sroundup_lwork(lwmin)
642*
643 RETURN
644*
645* End of CHBEVX_2STAGE
646*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
subroutine cgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
CGEMV
Definition cgemv.f:160
subroutine chetrd_hb2st(stage1, vect, uplo, n, kd, ab, ldab, d, e, hous, lhous, work, lwork, info)
CHETRD_HB2ST reduces a complex Hermitian band matrix A to real symmetric tridiagonal form T
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 clanhb(norm, uplo, n, k, ab, ldab, work)
CLANHB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition clanhb.f:132
subroutine clascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition clascl.f:143
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
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
Here is the call graph for this function:
Here is the caller graph for this function: