LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ ssyevx_2stage()

subroutine ssyevx_2stage ( character  JOBZ,
character  RANGE,
character  UPLO,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
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 
)

SSYEVX_2STAGE computes the eigenvalues and, optionally, the left and/or right eigenvectors for SY matrices

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

Purpose:
 SSYEVX_2STAGE computes selected eigenvalues and, optionally, eigenvectors
 of a real symmetric 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 REAL array, dimension (LDA, N)
          On entry, the symmetric 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 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 (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 + 3*N
                                             = N*KD + N*max(KD+1,FACTOPTNB) 
                                               + max(2*KD*KD, KD*NTHREADS) 
                                               + (KD+1)*N + 3*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]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 302 of file ssyevx_2stage.f.

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