LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ 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.
Date
June 2016
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 329 of file chbevx_2stage.f.

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