LAPACK  3.7.0
LAPACK: Linear Algebra PACKage
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.7.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 ilaenv
340  REAL slamch, slansy
341  EXTERNAL lsame, ilaenv, slamch, slansy
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 = ilaenv( 17, 'SSYTRD_2STAGE', jobz, n, -1, -1, -1 )
397  ib = ilaenv( 18, 'SSYTRD_2STAGE', jobz, n, kd, -1, -1 )
398  lhtrd = ilaenv( 19, 'SSYTRD_2STAGE', jobz, n, kd, ib, -1 )
399  lwtrd = ilaenv( 20, 'SSYTRD_2STAGE', jobz, n, kd, ib, -1 )
400  lwmin = max( 8*n, 3*n + lhtrd + lwtrd )
401  work( 1 ) = lwmin
402  END IF
403 *
404  IF( lwork.LT.lwmin .AND. .NOT.lquery )
405  $ info = -17
406  END IF
407 *
408  IF( info.NE.0 ) THEN
409  CALL xerbla( 'SSYEVX_2STAGE', -info )
410  RETURN
411  ELSE IF( lquery ) THEN
412  RETURN
413  END IF
414 *
415 * Quick return if possible
416 *
417  m = 0
418  IF( n.EQ.0 ) THEN
419  RETURN
420  END IF
421 *
422  IF( n.EQ.1 ) THEN
423  IF( alleig .OR. indeig ) THEN
424  m = 1
425  w( 1 ) = a( 1, 1 )
426  ELSE
427  IF( vl.LT.a( 1, 1 ) .AND. vu.GE.a( 1, 1 ) ) THEN
428  m = 1
429  w( 1 ) = a( 1, 1 )
430  END IF
431  END IF
432  IF( wantz )
433  $ z( 1, 1 ) = one
434  RETURN
435  END IF
436 *
437 * Get machine constants.
438 *
439  safmin = slamch( 'Safe minimum' )
440  eps = slamch( 'Precision' )
441  smlnum = safmin / eps
442  bignum = one / smlnum
443  rmin = sqrt( smlnum )
444  rmax = min( sqrt( bignum ), one / sqrt( sqrt( safmin ) ) )
445 *
446 * Scale matrix to allowable range, if necessary.
447 *
448  iscale = 0
449  abstll = abstol
450  IF( valeig ) THEN
451  vll = vl
452  vuu = vu
453  END IF
454  anrm = slansy( 'M', uplo, n, a, lda, work )
455  IF( anrm.GT.zero .AND. anrm.LT.rmin ) THEN
456  iscale = 1
457  sigma = rmin / anrm
458  ELSE IF( anrm.GT.rmax ) THEN
459  iscale = 1
460  sigma = rmax / anrm
461  END IF
462  IF( iscale.EQ.1 ) THEN
463  IF( lower ) THEN
464  DO 10 j = 1, n
465  CALL sscal( n-j+1, sigma, a( j, j ), 1 )
466  10 CONTINUE
467  ELSE
468  DO 20 j = 1, n
469  CALL sscal( j, sigma, a( 1, j ), 1 )
470  20 CONTINUE
471  END IF
472  IF( abstol.GT.0 )
473  $ abstll = abstol*sigma
474  IF( valeig ) THEN
475  vll = vl*sigma
476  vuu = vu*sigma
477  END IF
478  END IF
479 *
480 * Call SSYTRD_2STAGE to reduce symmetric matrix to tridiagonal form.
481 *
482  indtau = 1
483  inde = indtau + n
484  indd = inde + n
485  indhous = indd + n
486  indwrk = indhous + lhtrd
487  llwork = lwork - indwrk + 1
488 *
489  CALL ssytrd_2stage( jobz, uplo, n, a, lda, work( indd ),
490  $ work( inde ), work( indtau ), work( indhous ),
491  $ lhtrd, work( indwrk ), llwork, iinfo )
492 *
493 * If all eigenvalues are desired and ABSTOL is less than or equal to
494 * zero, then call SSTERF or SORGTR and SSTEQR. If this fails for
495 * some eigenvalue, then try SSTEBZ.
496 *
497  test = .false.
498  IF( indeig ) THEN
499  IF( il.EQ.1 .AND. iu.EQ.n ) THEN
500  test = .true.
501  END IF
502  END IF
503  IF( ( alleig .OR. test ) .AND. ( abstol.LE.zero ) ) THEN
504  CALL scopy( n, work( indd ), 1, w, 1 )
505  indee = indwrk + 2*n
506  IF( .NOT.wantz ) THEN
507  CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
508  CALL ssterf( n, w, work( indee ), info )
509  ELSE
510  CALL slacpy( 'A', n, n, a, lda, z, ldz )
511  CALL sorgtr( uplo, n, z, ldz, work( indtau ),
512  $ work( indwrk ), llwork, iinfo )
513  CALL scopy( n-1, work( inde ), 1, work( indee ), 1 )
514  CALL ssteqr( jobz, n, w, work( indee ), z, ldz,
515  $ work( indwrk ), info )
516  IF( info.EQ.0 ) THEN
517  DO 30 i = 1, n
518  ifail( i ) = 0
519  30 CONTINUE
520  END IF
521  END IF
522  IF( info.EQ.0 ) THEN
523  m = n
524  GO TO 40
525  END IF
526  info = 0
527  END IF
528 *
529 * Otherwise, call SSTEBZ and, if eigenvectors are desired, SSTEIN.
530 *
531  IF( wantz ) THEN
532  order = 'B'
533  ELSE
534  order = 'E'
535  END IF
536  indibl = 1
537  indisp = indibl + n
538  indiwo = indisp + n
539  CALL sstebz( range, order, n, vll, vuu, il, iu, abstll,
540  $ work( indd ), work( inde ), m, nsplit, w,
541  $ iwork( indibl ), iwork( indisp ), work( indwrk ),
542  $ iwork( indiwo ), info )
543 *
544  IF( wantz ) THEN
545  CALL sstein( n, work( indd ), work( inde ), m, w,
546  $ iwork( indibl ), iwork( indisp ), z, ldz,
547  $ work( indwrk ), iwork( indiwo ), ifail, info )
548 *
549 * Apply orthogonal matrix used in reduction to tridiagonal
550 * form to eigenvectors returned by SSTEIN.
551 *
552  indwkn = inde
553  llwrkn = lwork - indwkn + 1
554  CALL sormtr( 'L', uplo, 'N', n, m, a, lda, work( indtau ), z,
555  $ ldz, work( indwkn ), llwrkn, iinfo )
556  END IF
557 *
558 * If matrix was scaled, then rescale eigenvalues appropriately.
559 *
560  40 CONTINUE
561  IF( iscale.EQ.1 ) THEN
562  IF( info.EQ.0 ) THEN
563  imax = m
564  ELSE
565  imax = info - 1
566  END IF
567  CALL sscal( imax, one / sigma, w, 1 )
568  END IF
569 *
570 * If eigenvalues are not in order, then sort them, along with
571 * eigenvectors.
572 *
573  IF( wantz ) THEN
574  DO 60 j = 1, m - 1
575  i = 0
576  tmp1 = w( j )
577  DO 50 jj = j + 1, m
578  IF( w( jj ).LT.tmp1 ) THEN
579  i = jj
580  tmp1 = w( jj )
581  END IF
582  50 CONTINUE
583 *
584  IF( i.NE.0 ) THEN
585  itmp1 = iwork( indibl+i-1 )
586  w( i ) = w( j )
587  iwork( indibl+i-1 ) = iwork( indibl+j-1 )
588  w( j ) = tmp1
589  iwork( indibl+j-1 ) = itmp1
590  CALL sswap( n, z( 1, i ), 1, z( 1, j ), 1 )
591  IF( info.NE.0 ) THEN
592  itmp1 = ifail( i )
593  ifail( i ) = ifail( j )
594  ifail( j ) = itmp1
595  END IF
596  END IF
597  60 CONTINUE
598  END IF
599 *
600 * Set WORK(1) to optimal workspace size.
601 *
602  work( 1 ) = lwmin
603 *
604  RETURN
605 *
606 * End of SSYEVX_2STAGE
607 *
subroutine sormtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMTR
Definition: sormtr.f:174
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:133
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
subroutine sstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
SSTEIN
Definition: sstein.f:176
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
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
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:53
subroutine ssytrd_2stage(VECT, UPLO, N, A, LDA, D, E, TAU, HOUS2, LHOUS2, WORK, LWORK, INFO)
SSYTRD_2STAGE
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
subroutine ssterf(N, D, E, INFO)
SSTERF
Definition: ssterf.f:88
subroutine sorgtr(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO)
SORGTR
Definition: sorgtr.f:125
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
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, or the element of largest absolute value of a real symmetric matrix.
Definition: slansy.f:124

Here is the call graph for this function:

Here is the caller graph for this function: