LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ ssbevx_2stage()

subroutine ssbevx_2stage ( character  JOBZ,
character  RANGE,
character  UPLO,
integer  N,
integer  KD,
real, dimension( ldab, * )  AB,
integer  LDAB,
real, dimension( ldq, * )  Q,
integer  LDQ,
real  VL,
real  VU,
integer  IL,
integer  IU,
real  ABSTOL,
integer  M,
real, dimension( * )  W,
real, dimension( ldz, * )  Z,
integer  LDZ,
real, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
integer, dimension( * )  IFAIL,
integer  INFO 
)

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

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

Purpose:
 SSBEVX_2STAGE computes selected eigenvalues and, optionally, eigenvectors
 of a real symmetric 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 REAL array, dimension (LDAB, N)
          On entry, the upper or lower triangle of the symmetric 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.  If UPLO = 'U', the first
          superdiagonal and the diagonal of the tridiagonal matrix T
          are returned in rows KD and KD+1 of AB, and if UPLO = 'L',
          the diagonal and first subdiagonal of T are returned in the
          first two rows of AB.
[in]LDAB
          LDAB is INTEGER
          The leading dimension of the array AB.  LDAB >= KD + 1.
[out]Q
          Q is REAL array, dimension (LDQ, N)
          If JOBZ = 'V', the N-by-N orthogonal 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 REAL 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 REAL 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, 7*N, dimension) where
                                   dimension = (2KD+1)*N + KD*NTHREADS + 2*N
                                   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 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]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 324 of file ssbevx_2stage.f.

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