LAPACK  3.10.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.
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 319 of file ssbevx_2stage.f.

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