 LAPACK  3.8.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

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.```
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).
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 *
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 ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:133
subroutine ssytrd_2stage(VECT, UPLO, N, A, LDA, D, E, TAU, HOUS2, LHOUS2, WORK, LWORK, INFO)
SSYTRD_2STAGE
subroutine sormtr(SIDE, UPLO, TRANS, M, N, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMTR
Definition: sormtr.f:174
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 sstein(N, D, E, M, W, IBLOCK, ISPLIT, Z, LDZ, WORK, IWORK, IFAIL, INFO)
SSTEIN
Definition: sstein.f:176
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:81
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:84
integer function ilaenv2stage(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV2STAGE
Definition: ilaenv2stage.f:151
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:84
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: